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.
iceitd.F90 in NEMO/releases/r4.0/r4.0-HEAD/src/ICE – NEMO

source: NEMO/releases/r4.0/r4.0-HEAD/src/ICE/iceitd.F90

Last change on this file was 15045, checked in by clem, 3 years ago

4.0-HEAD: bug found by ecmwf: do not do ice remapping if, by some mystery and in very rare occasion, hR=hL=hice.
  • Property svn:keywords set to Id
File size: 38.1 KB
Line 
1MODULE iceitd
2   !!======================================================================
3   !!                       ***  MODULE iceitd ***
4   !!   sea-ice : ice thickness distribution
5   !!======================================================================
6   !! History :  3.0  !  2005-12  (M. Vancoppenolle) original code (based on CICE)
7   !!            4.0  !  2018     (many people)      SI3 [aka Sea Ice cube]
8   !!----------------------------------------------------------------------
9#if defined key_si3
10   !!----------------------------------------------------------------------
11   !!   'key_si3'                                       SI3 sea-ice model
12   !!----------------------------------------------------------------------
13   !!   ice_itd_rem   : redistribute ice thicknesses after thermo growth and melt
14   !!   itd_glinear   : build g(h) satisfying area and volume constraints
15   !!   itd_shiftice  : shift ice across category boundaries, conserving everything
16   !!   ice_itd_reb   : rebin ice thicknesses into bounded categories
17   !!   ice_itd_init  : read ice thicknesses mean and min from namelist
18   !!----------------------------------------------------------------------
19   USE dom_oce        ! ocean domain
20   USE phycst         ! physical constants
21   USE ice1D          ! sea-ice: thermodynamic variables
22   USE ice            ! sea-ice: variables
23   USE icevar         ! sea-ice: operations
24   USE icectl         ! sea-ice: conservation tests
25   USE icetab         ! sea-ice: convert 1D<=>2D
26   !
27   USE in_out_manager ! I/O manager
28   USE lib_mpp        ! MPP library
29   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
30   USE prtctl         ! Print control
31   USE timing         ! Timing
32
33   IMPLICIT NONE
34   PRIVATE
35
36   PUBLIC   ice_itd_init  ! called in icestp
37   PUBLIC   ice_itd_rem   ! called in icethd
38   PUBLIC   ice_itd_reb   ! called in icecor
39
40   INTEGER            ::   nice_catbnd     ! choice of the type of ice category function
41   !                                       ! associated indices:
42   INTEGER, PARAMETER ::   np_cathfn = 1   ! categories defined by a function
43   INTEGER, PARAMETER ::   np_catusr = 2   ! categories defined by the user
44   !
45   !                                           !! ** namelist (namitd) **
46   LOGICAL                    ::   ln_cat_hfn   ! ice categories are defined by function like rn_himean**(-0.05)
47   REAL(wp)                   ::   rn_himean    ! mean thickness of the domain
48   LOGICAL                    ::   ln_cat_usr   ! ice categories are defined by rn_catbnd
49   REAL(wp), DIMENSION(0:100) ::   rn_catbnd    ! ice categories bounds
50   REAL(wp)                   ::   rn_himax     ! maximum ice thickness allowed
51   !
52   !!----------------------------------------------------------------------
53   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
54   !! $Id$
55   !! Software governed by the CeCILL license (see ./LICENSE)
56   !!----------------------------------------------------------------------
57CONTAINS
58
59   SUBROUTINE ice_itd_rem( kt )
60      !!------------------------------------------------------------------
61      !!                ***  ROUTINE ice_itd_rem ***
62      !!
63      !! ** Purpose :   computes the redistribution of ice thickness
64      !!                after thermodynamic growth of ice thickness
65      !!
66      !! ** Method  :   Linear remapping
67      !!
68      !! References :   W.H. Lipscomb, JGR 2001
69      !!------------------------------------------------------------------
70      INTEGER , INTENT (in) ::   kt      ! Ocean time step
71      !
72      INTEGER  ::   ji, jj, jl, jcat     ! dummy loop index
73      INTEGER  ::   ipti                 ! local integer
74      REAL(wp) ::   zx1, zwk1, zdh0, zetamin, zdamax   ! local scalars
75      REAL(wp) ::   zx2, zwk2, zda0, zetamax           !   -      -
76      REAL(wp) ::   zx3       
77      REAL(wp) ::   zslope          ! used to compute local thermodynamic "speeds"
78      !
79      INTEGER , DIMENSION(jpij)       ::   iptidx          ! compute remapping or not
80      INTEGER , DIMENSION(jpij,jpl-1) ::   jdonor          ! donor category index
81      REAL(wp), DIMENSION(jpij,jpl)   ::   zdhice          ! ice thickness increment
82      REAL(wp), DIMENSION(jpij,jpl)   ::   g0, g1          ! coefficients for fitting the line of the ITD
83      REAL(wp), DIMENSION(jpij,jpl)   ::   hL, hR          ! left and right boundary for the ITD for each thickness
84      REAL(wp), DIMENSION(jpij,jpl-1) ::   zdaice, zdvice  ! local increment of ice area and volume
85      REAL(wp), DIMENSION(jpij)       ::   zhb0, zhb1      ! category boundaries for thinnes categories
86      REAL(wp), DIMENSION(jpij,0:jpl) ::   zhbnew          ! new boundaries of ice categories
87      !!------------------------------------------------------------------
88      IF( ln_timing )   CALL timing_start('iceitd_rem')
89
90      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_itd_rem: remapping ice thickness distribution' 
91
92      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'iceitd_rem', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
93      IF( ln_icediachk )   CALL ice_cons2D  (0, 'iceitd_rem',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft)
94
95      !-----------------------------------------------------------------------------------------------
96      !  1) Identify grid cells with ice
97      !-----------------------------------------------------------------------------------------------
98      at_i(:,:) = SUM( a_i, dim=3 )
99      !
100      npti = 0   ;   nptidx(:) = 0
101      DO jj = 1, jpj
102         DO ji = 1, jpi
103            IF ( at_i(ji,jj) > epsi10 ) THEN
104               npti = npti + 1
105               nptidx( npti ) = (jj - 1) * jpi + ji
106            ENDIF
107         END DO
108      END DO
109     
110      !-----------------------------------------------------------------------------------------------
111      !  2) Compute new category boundaries
112      !-----------------------------------------------------------------------------------------------
113      IF( npti > 0 ) THEN
114         !
115         zdhice(:,:) = 0._wp
116         zhbnew(:,:) = 0._wp
117         !
118         CALL tab_3d_2d( npti, nptidx(1:npti), h_i_2d (1:npti,1:jpl), h_i   )
119         CALL tab_3d_2d( npti, nptidx(1:npti), h_ib_2d(1:npti,1:jpl), h_i_b )
120         CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i   )
121         CALL tab_3d_2d( npti, nptidx(1:npti), a_ib_2d(1:npti,1:jpl), a_i_b )
122         !
123         DO jl = 1, jpl
124            ! Compute thickness change in each ice category
125            DO ji = 1, npti
126               IF( a_i_2d(ji,jl) > epsi10 )   zdhice(ji,jl) = h_i_2d(ji,jl) - h_ib_2d(ji,jl)
127            END DO
128         END DO
129         !
130         ! --- New boundaries for category 1:jpl-1 --- !
131         DO jl = 1, jpl - 1
132            !
133            DO ji = 1, npti
134               !
135               ! --- New boundary: Hn* = Hn + Fn*dt --- !
136               !     Fn*dt = ( fn + (fn+1 - fn)/(hn+1 - hn) * (Hn - hn) ) * dt = zdhice + zslope * (Hmax - h_i_b)
137               !
138               IF    ( a_ib_2d(ji,jl) >  epsi10 .AND. a_ib_2d(ji,jl+1) >  epsi10 ) THEN   ! a(jl+1) & a(jl) /= 0
139                  zslope        = ( zdhice(ji,jl+1) - zdhice(ji,jl) ) / ( h_ib_2d(ji,jl+1) - h_ib_2d(ji,jl) )
140                  zhbnew(ji,jl) = hi_max(jl) + zdhice(ji,jl) + zslope * ( hi_max(jl) - h_ib_2d(ji,jl) )
141               ELSEIF( a_ib_2d(ji,jl) >  epsi10 .AND. a_ib_2d(ji,jl+1) <= epsi10 ) THEN   ! a(jl+1)=0 => Hn* = Hn + fn*dt
142                  zhbnew(ji,jl) = hi_max(jl) + zdhice(ji,jl)
143               ELSEIF( a_ib_2d(ji,jl) <= epsi10 .AND. a_ib_2d(ji,jl+1) >  epsi10 ) THEN   ! a(jl)=0 => Hn* = Hn + fn+1*dt
144                  zhbnew(ji,jl) = hi_max(jl) + zdhice(ji,jl+1)
145               ELSE                                                                       ! a(jl+1) & a(jl) = 0
146                  zhbnew(ji,jl) = hi_max(jl)
147               ENDIF
148               !
149               ! --- 2 conditions for remapping --- !
150               ! 1) hn(t+1)+espi < Hn* < hn+1(t+1)-epsi               
151               !    Note: hn(t+1) must not be too close to either HR or HL otherwise a division by nearly 0 is possible
152               !          in itd_glinear in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice)
153               IF( a_i_2d(ji,jl  ) > epsi10 .AND. h_i_2d(ji,jl  ) > ( zhbnew(ji,jl) - epsi10 ) )   nptidx(ji) = 0
154               IF( a_i_2d(ji,jl+1) > epsi10 .AND. h_i_2d(ji,jl+1) < ( zhbnew(ji,jl) + epsi10 ) )   nptidx(ji) = 0
155               !
156               ! 2) Hn-1 < Hn* < Hn+1 
157               IF( zhbnew(ji,jl) < hi_max(jl-1) )   nptidx(ji) = 0
158               IF( zhbnew(ji,jl) > hi_max(jl+1) )   nptidx(ji) = 0
159               !
160            END DO
161         END DO
162         !
163         ! --- New boundaries for category jpl --- !
164         DO ji = 1, npti
165            IF( a_i_2d(ji,jpl) > epsi10 ) THEN
166               zhbnew(ji,jpl) = MAX( hi_max(jpl-1), 3._wp * h_i_2d(ji,jpl) - 2._wp * zhbnew(ji,jpl-1) )
167            ELSE
168               zhbnew(ji,jpl) = hi_max(jpl) 
169            ENDIF
170            !
171            ! --- 1 additional condition for remapping (1st category) --- !
172            ! H0+epsi < h1(t) < H1-epsi
173            !    h1(t) must not be too close to either HR or HL otherwise a division by nearly 0 is possible
174            !    in itd_glinear in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice)
175            IF( h_ib_2d(ji,1) < ( hi_max(0) + epsi10 ) )   nptidx(ji) = 0
176            IF( h_ib_2d(ji,1) > ( hi_max(1) - epsi10 ) )   nptidx(ji) = 0
177         END DO
178         !
179         !-----------------------------------------------------------------------------------------------
180         !  3) Identify cells where remapping
181         !-----------------------------------------------------------------------------------------------
182         ipti = 0   ;   iptidx(:) = 0
183         DO ji = 1, npti
184            IF( nptidx(ji) /= 0 ) THEN
185               ipti = ipti + 1
186               iptidx(ipti)   = nptidx(ji)
187               zhbnew(ipti,:) = zhbnew(ji,:)  ! adjust zhbnew to new indices
188            ENDIF
189         END DO
190         nptidx(:) = iptidx(:)
191         npti      = ipti
192         !
193      ENDIF
194   
195      !-----------------------------------------------------------------------------------------------
196      !  4) Compute g(h)
197      !-----------------------------------------------------------------------------------------------
198      IF( npti > 0 ) THEN
199         !
200         zhb0(:) = hi_max(0)   ;   zhb1(:) = hi_max(1)
201         g0(:,:) = 0._wp       ;   g1(:,:) = 0._wp 
202         hL(:,:) = 0._wp       ;   hR(:,:) = 0._wp 
203         !
204         DO jl = 1, jpl
205            !
206            CALL tab_2d_1d( npti, nptidx(1:npti), h_ib_1d(1:npti), h_i_b(:,:,jl) )
207            CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i  (:,:,jl) )
208            CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i  (:,:,jl) )
209            CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d (1:npti), v_i  (:,:,jl) )
210            !
211            IF( jl == 1 ) THEN
212               ! 
213               ! --- g(h) for category 1 --- !
214               CALL itd_glinear( zhb0(1:npti)  , zhb1(1:npti)  , h_ib_1d(1:npti)  , a_i_1d(1:npti)  ,  &   ! in
215                  &              g0  (1:npti,1), g1  (1:npti,1), hL     (1:npti,1), hR    (1:npti,1)   )   ! out
216               !
217               ! Area lost due to melting of thin ice
218               DO ji = 1, npti
219                  !
220                  IF( a_i_1d(ji) > epsi10 ) THEN
221                     !
222                     zdh0 =  h_i_1d(ji) - h_ib_1d(ji)               
223                     IF( zdh0 < 0.0 ) THEN      ! remove area from category 1
224                        zdh0 = MIN( -zdh0, hi_max(1) )
225                        !Integrate g(1) from 0 to dh0 to estimate area melted
226                        zetamax = MIN( zdh0, hR(ji,1) ) - hL(ji,1)
227                        !
228                        IF( zetamax > 0.0 ) THEN
229                           zx1    = zetamax
230                           zx2    = 0.5 * zetamax * zetamax 
231                           zda0   = g1(ji,1) * zx2 + g0(ji,1) * zx1                ! ice area removed
232                           zdamax = a_i_1d(ji) * (1.0 - h_i_1d(ji) / h_ib_1d(ji) ) ! Constrain new thickness <= h_i               
233                           zda0   = MIN( zda0, zdamax )                            ! ice area lost due to melting of thin ice (zdamax > 0)
234                           ! Remove area, conserving volume
235                           h_i_1d(ji) = h_i_1d(ji) * a_i_1d(ji) / ( a_i_1d(ji) - zda0 )
236                           a_i_1d(ji) = a_i_1d(ji) - zda0
237                           v_i_1d(ji) = a_i_1d(ji) * h_i_1d(ji) ! useless ?
238                        ENDIF
239                        !
240                     ELSE ! if ice accretion zdh0 > 0
241                        ! zhbnew was 0, and is shifted to the right to account for thin ice growth in openwater (F0 = f1)
242                        zhbnew(ji,0) = MIN( zdh0, hi_max(1) ) 
243                     ENDIF
244                     !
245                  ENDIF
246                  !
247               END DO
248               !
249               CALL tab_1d_2d( npti, nptidx(1:npti), h_i_1d(1:npti), h_i(:,:,jl) )
250               CALL tab_1d_2d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl) )
251               CALL tab_1d_2d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl) )
252               !
253            ENDIF ! jl=1
254            !
255            ! --- g(h) for each thickness category --- ! 
256            CALL itd_glinear( zhbnew(1:npti,jl-1), zhbnew(1:npti,jl), h_i_1d(1:npti)   , a_i_1d(1:npti)   ,  &   ! in
257               &              g0    (1:npti,jl  ), g1    (1:npti,jl), hL    (1:npti,jl), hR    (1:npti,jl)   )   ! out
258            !
259         END DO
260         
261         !-----------------------------------------------------------------------------------------------
262         !  5) Compute area and volume to be shifted across each boundary (Eq. 18)
263         !-----------------------------------------------------------------------------------------------
264         DO jl = 1, jpl - 1
265            !
266            DO ji = 1, npti
267               !
268               ! left and right integration limits in eta space
269               IF (zhbnew(ji,jl) > hi_max(jl)) THEN ! Hn* > Hn => transfer from jl to jl+1
270                  zetamin = MAX( hi_max(jl)   , hL(ji,jl) ) - hL(ji,jl)   ! hi_max(jl) - hL
271                  zetamax = MIN( zhbnew(ji,jl), hR(ji,jl) ) - hL(ji,jl)   ! hR - hL
272                  jdonor(ji,jl) = jl
273               ELSE                                 ! Hn* <= Hn => transfer from jl+1 to jl
274                  zetamin = 0.0
275                  zetamax = MIN( hi_max(jl), hR(ji,jl+1) ) - hL(ji,jl+1)  ! hi_max(jl) - hL
276                  jdonor(ji,jl) = jl + 1
277               ENDIF
278               zetamax = MAX( zetamax, zetamin ) ! no transfer if etamax < etamin
279               !
280               zx1  = zetamax - zetamin
281               zwk1 = zetamin * zetamin
282               zwk2 = zetamax * zetamax
283               zx2  = 0.5 * ( zwk2 - zwk1 )
284               zwk1 = zwk1 * zetamin
285               zwk2 = zwk2 * zetamax
286               zx3  = 1.0 / 3.0 * ( zwk2 - zwk1 )
287               jcat = jdonor(ji,jl)
288               zdaice(ji,jl) = g1(ji,jcat)*zx2 + g0(ji,jcat)*zx1
289               zdvice(ji,jl) = g1(ji,jcat)*zx3 + g0(ji,jcat)*zx2 + zdaice(ji,jl)*hL(ji,jcat)
290               !
291            END DO
292         END DO
293         
294         !----------------------------------------------------------------------------------------------
295         ! 6) Shift ice between categories
296         !----------------------------------------------------------------------------------------------
297         CALL itd_shiftice ( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) )
298         
299         !----------------------------------------------------------------------------------------------
300         ! 7) Make sure h_i >= minimum ice thickness hi_min
301         !----------------------------------------------------------------------------------------------
302         CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i (:,:,1) )
303         CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i (:,:,1) )
304         CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d(1:npti), a_ip(:,:,1) )
305         !
306         DO ji = 1, npti
307            IF ( a_i_1d(ji) > epsi10 .AND. h_i_1d(ji) < rn_himin ) THEN
308               a_i_1d(ji) = a_i_1d(ji) * h_i_1d(ji) / rn_himin 
309               IF( ln_pnd_LEV )   a_ip_1d(ji) = a_ip_1d(ji) * h_i_1d(ji) / rn_himin
310               h_i_1d(ji) = rn_himin
311            ENDIF
312         END DO
313         !
314         CALL tab_1d_2d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i (:,:,1) )
315         CALL tab_1d_2d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i (:,:,1) )
316         CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d(1:npti), a_ip(:,:,1) )
317         !
318      ENDIF
319      !
320      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'iceitd_rem', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
321      IF( ln_icediachk )   CALL ice_cons2D  (1, 'iceitd_rem',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft)
322      IF( ln_timing    )   CALL timing_stop ('iceitd_rem')
323      !
324   END SUBROUTINE ice_itd_rem
325
326
327   SUBROUTINE itd_glinear( HbL, Hbr, phice, paice, pg0, pg1, phL, phR )
328      !!------------------------------------------------------------------
329      !!                ***  ROUTINE itd_glinear ***
330      !!
331      !! ** Purpose :   build g(h) satisfying area and volume constraints (Eq. 6 and 9)
332      !!
333      !! ** Method  :   g(h) is linear and written as: g(eta) = g1(eta) + g0
334      !!                with eta = h - HL
335      !!------------------------------------------------------------------
336      REAL(wp), DIMENSION(:), INTENT(in   ) ::   HbL, HbR      ! left and right category boundaries
337      REAL(wp), DIMENSION(:), INTENT(in   ) ::   phice, paice  ! ice thickness and concentration
338      REAL(wp), DIMENSION(:), INTENT(inout) ::   pg0, pg1      ! coefficients in linear equation for g(eta)
339      REAL(wp), DIMENSION(:), INTENT(inout) ::   phL, phR      ! min and max value of range over which g(h) > 0
340      !
341      INTEGER  ::   ji           ! horizontal indices
342      REAL(wp) ::   z1_3 , z2_3  ! 1/3 , 2/3
343      REAL(wp) ::   zh13         ! HbL + 1/3 * (HbR - HbL)
344      REAL(wp) ::   zh23         ! HbL + 2/3 * (HbR - HbL)
345      REAL(wp) ::   zdhr         ! 1 / (hR - hL)
346      REAL(wp) ::   zwk1, zwk2   ! temporary variables
347      !!------------------------------------------------------------------
348      !
349      z1_3 = 1._wp / 3._wp
350      z2_3 = 2._wp / 3._wp
351      !
352      DO ji = 1, npti
353         !
354         IF( paice(ji) > epsi10  .AND. phice(ji) > epsi10 )  THEN
355            !
356            ! Initialize hL and hR
357            phL(ji) = HbL(ji)
358            phR(ji) = HbR(ji)
359            !
360            ! Change hL or hR if hice falls outside central third of range,
361            ! so that hice is in the central third of the range [HL HR]
362            zh13 = z1_3 * ( 2._wp * phL(ji) +         phR(ji) )
363            zh23 = z1_3 * (         phL(ji) + 2._wp * phR(ji) )
364            !
365            IF    ( phice(ji) < zh13 ) THEN   ;   phR(ji) = 3._wp * phice(ji) - 2._wp * phL(ji) ! move HR to the left
366            ELSEIF( phice(ji) > zh23 ) THEN   ;   phL(ji) = 3._wp * phice(ji) - 2._wp * phR(ji) ! move HL to the right
367            ENDIF
368            !
369            ! Compute coefficients of g(eta) = g0 + g1*eta
370            IF( phR(ji) > phL(ji) ) THEN   ;   zdhr = 1._wp / (phR(ji) - phL(ji))
371            ELSE                           ;   zdhr = 0._wp ! if hR=hL=hice => no remapping
372            ENDIF
373            !!zdhr = 1._wp / (phR(ji) - phL(ji))
374            zwk1 = 6._wp * paice(ji) * zdhr
375            zwk2 = ( phice(ji) - phL(ji) ) * zdhr
376            pg0(ji) = zwk1 * ( z2_3 - zwk2 )                    ! Eq. 14
377            pg1(ji) = 2._wp * zdhr * zwk1 * ( zwk2 - 0.5_wp )   ! Eq. 14
378            !
379         ELSE  ! remap_flag = .false. or a_i < epsi10
380            phL(ji) = 0._wp
381            phR(ji) = 0._wp
382            pg0(ji) = 0._wp
383            pg1(ji) = 0._wp
384         ENDIF
385         !
386      END DO
387      !
388   END SUBROUTINE itd_glinear
389
390
391   SUBROUTINE itd_shiftice( kdonor, pdaice, pdvice )
392      !!------------------------------------------------------------------
393      !!                ***  ROUTINE itd_shiftice ***
394      !!
395      !! ** Purpose :   shift ice across category boundaries, conserving everything
396      !!              ( area, volume, energy, age*vol, and mass of salt )
397      !!------------------------------------------------------------------
398      INTEGER , DIMENSION(:,:), INTENT(in) ::   kdonor   ! donor category index
399      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pdaice   ! ice area transferred across boundary
400      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pdvice   ! ice volume transferred across boundary
401      !
402      INTEGER  ::   ji, jl, jk         ! dummy loop indices
403      INTEGER  ::   jl2, jl1           ! local integers
404      REAL(wp) ::   ztrans             ! ice/snow transferred
405      REAL(wp), DIMENSION(jpij)            ::   zworka, zworkv   ! workspace
406      REAL(wp), DIMENSION(jpij,jpl)        ::   zaTsfn           !  -    -
407      REAL(wp), DIMENSION(jpij,nlay_i,jpl) ::   ze_i_2d
408      REAL(wp), DIMENSION(jpij,nlay_s,jpl) ::   ze_s_2d
409      !!------------------------------------------------------------------
410         
411      CALL tab_3d_2d( npti, nptidx(1:npti), h_i_2d (1:npti,1:jpl), h_i  )
412      CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i  )
413      CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i  )
414      CALL tab_3d_2d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s  )
415      CALL tab_3d_2d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i )
416      CALL tab_3d_2d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i )
417      CALL tab_3d_2d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip )
418      CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip )
419      CALL tab_3d_2d( npti, nptidx(1:npti), v_il_2d(1:npti,1:jpl), v_il )
420      CALL tab_3d_2d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su )
421      DO jl = 1, jpl
422         DO jk = 1, nlay_s
423            CALL tab_2d_1d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) )
424         END DO
425         DO jk = 1, nlay_i
426            CALL tab_2d_1d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
427         END DO
428      END DO
429      ! to correct roundoff errors on a_i
430      CALL tab_2d_1d( npti, nptidx(1:npti), rn_amax_1d(1:npti), rn_amax_2d )
431
432      !----------------------------------------------------------------------------------------------
433      ! 1) Define a variable equal to a_i*T_su
434      !----------------------------------------------------------------------------------------------
435      DO jl = 1, jpl
436         DO ji = 1, npti
437            zaTsfn(ji,jl) = a_i_2d(ji,jl) * t_su_2d(ji,jl)
438         END DO
439      END DO
440     
441      !-------------------------------------------------------------------------------
442      ! 2) Transfer volume and energy between categories
443      !-------------------------------------------------------------------------------
444      DO jl = 1, jpl - 1
445         DO ji = 1, npti
446            !
447            jl1 = kdonor(ji,jl)
448            !
449            IF( jl1 > 0 ) THEN
450               !
451               IF ( jl1 == jl  ) THEN   ;   jl2 = jl1+1
452               ELSE                     ;   jl2 = jl 
453               ENDIF
454               !
455               IF( v_i_2d(ji,jl1) >= epsi10 ) THEN   ;   zworkv(ji) = pdvice(ji,jl) / v_i_2d(ji,jl1)
456               ELSE                                  ;   zworkv(ji) = 0._wp
457               ENDIF
458               IF( a_i_2d(ji,jl1) >= epsi10 ) THEN   ;   zworka(ji) = pdaice(ji,jl) / a_i_2d(ji,jl1)
459               ELSE                                  ;   zworka(ji) = 0._wp
460               ENDIF
461               !
462               a_i_2d(ji,jl1) = a_i_2d(ji,jl1) - pdaice(ji,jl)       ! Ice areas
463               a_i_2d(ji,jl2) = a_i_2d(ji,jl2) + pdaice(ji,jl)
464               !
465               v_i_2d(ji,jl1) = v_i_2d(ji,jl1) - pdvice(ji,jl)       ! Ice volumes
466               v_i_2d(ji,jl2) = v_i_2d(ji,jl2) + pdvice(ji,jl)
467               !
468               ztrans         = v_s_2d(ji,jl1) * zworkv(ji)          ! Snow volumes
469               v_s_2d(ji,jl1) = v_s_2d(ji,jl1) - ztrans
470               v_s_2d(ji,jl2) = v_s_2d(ji,jl2) + ztrans 
471               !
472               ztrans          = oa_i_2d(ji,jl1) * zworka(ji)        ! Ice age
473               oa_i_2d(ji,jl1) = oa_i_2d(ji,jl1) - ztrans
474               oa_i_2d(ji,jl2) = oa_i_2d(ji,jl2) + ztrans
475               !
476               ztrans          = sv_i_2d(ji,jl1) * zworkv(ji)        ! Ice salinity
477               sv_i_2d(ji,jl1) = sv_i_2d(ji,jl1) - ztrans
478               sv_i_2d(ji,jl2) = sv_i_2d(ji,jl2) + ztrans
479               !
480               ztrans          = zaTsfn(ji,jl1) * zworka(ji)         ! Surface temperature
481               zaTsfn(ji,jl1)  = zaTsfn(ji,jl1) - ztrans
482               zaTsfn(ji,jl2)  = zaTsfn(ji,jl2) + ztrans
483               ! 
484               IF ( ln_pnd_LEV ) THEN
485                  ztrans          = a_ip_2d(ji,jl1) * zworka(ji)     ! Pond fraction
486                  a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - ztrans
487                  a_ip_2d(ji,jl2) = a_ip_2d(ji,jl2) + ztrans
488                  !                                             
489!!$                  ztrans          = v_ip_2d(ji,jl1) * zworka(ji)     ! Pond volume (also proportional to da/a)
490                  ztrans          = v_ip_2d(ji,jl1) * zworkv(ji)     ! Pond volume
491                  v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - ztrans
492                  v_ip_2d(ji,jl2) = v_ip_2d(ji,jl2) + ztrans
493                  !
494                  IF ( ln_pnd_lids ) THEN                            ! Pond lid volume
495!!$                     ztrans          = v_il_2d(ji,jl1) * zworka(ji)
496                     ztrans          = v_il_2d(ji,jl1) * zworkv(ji)
497                     v_il_2d(ji,jl1) = v_il_2d(ji,jl1) - ztrans
498                     v_il_2d(ji,jl2) = v_il_2d(ji,jl2) + ztrans
499                  ENDIF
500               ENDIF
501               !
502            ENDIF   ! jl1 >0
503         END DO
504         !
505         DO jk = 1, nlay_s         !--- Snow heat content
506            DO ji = 1, npti
507               !
508               jl1 = kdonor(ji,jl)
509               !
510               IF( jl1 > 0 ) THEN
511                  IF(jl1 == jl) THEN  ;  jl2 = jl+1
512                  ELSE                ;  jl2 = jl
513                  ENDIF
514                  ztrans             = ze_s_2d(ji,jk,jl1) * zworkv(ji)
515                  ze_s_2d(ji,jk,jl1) = ze_s_2d(ji,jk,jl1) - ztrans
516                  ze_s_2d(ji,jk,jl2) = ze_s_2d(ji,jk,jl2) + ztrans
517               ENDIF
518            END DO
519         END DO
520         !
521         DO jk = 1, nlay_i         !--- Ice heat content
522            DO ji = 1, npti
523               !
524               jl1 = kdonor(ji,jl)
525               !
526               IF( jl1 > 0 ) THEN
527                  IF(jl1 == jl) THEN  ;  jl2 = jl+1
528                  ELSE                ;  jl2 = jl
529                  ENDIF
530                  ztrans             = ze_i_2d(ji,jk,jl1) * zworkv(ji)
531                  ze_i_2d(ji,jk,jl1) = ze_i_2d(ji,jk,jl1) - ztrans
532                  ze_i_2d(ji,jk,jl2) = ze_i_2d(ji,jk,jl2) + ztrans
533               ENDIF
534            END DO
535         END DO
536         !
537      END DO                   ! boundaries, 1 to jpl-1
538
539      !-------------------
540      ! 3) roundoff errors
541      !-------------------
542      ! clem: The transfer between one category to another can lead to very small negative values (-1.e-20)
543      !       because of truncation error ( i.e. 1. - 1. /= 0 )
544      CALL ice_var_roundoff( a_i_2d, v_i_2d, v_s_2d, sv_i_2d, oa_i_2d, a_ip_2d, v_ip_2d, v_il_2d, ze_s_2d, ze_i_2d )
545
546      ! at_i must be <= rn_amax
547      zworka(1:npti) = SUM( a_i_2d(1:npti,:), dim=2 )
548      DO jl  = 1, jpl
549         WHERE( zworka(1:npti) > rn_amax_1d(1:npti) )   &
550            &   a_i_2d(1:npti,jl) = a_i_2d(1:npti,jl) * rn_amax_1d(1:npti) / zworka(1:npti)
551      END DO
552     
553      !-------------------------------------------------------------------------------
554      ! 4) Update ice thickness and temperature
555      !-------------------------------------------------------------------------------
556      WHERE( a_i_2d(1:npti,:) >= epsi20 )
557         h_i_2d (1:npti,:)  =  v_i_2d(1:npti,:) / a_i_2d(1:npti,:) 
558         t_su_2d(1:npti,:)  =  zaTsfn(1:npti,:) / a_i_2d(1:npti,:) 
559      ELSEWHERE
560         h_i_2d (1:npti,:)  = 0._wp
561         t_su_2d(1:npti,:)  = rt0
562      END WHERE
563      !
564      CALL tab_2d_3d( npti, nptidx(1:npti), h_i_2d (1:npti,1:jpl), h_i  )
565      CALL tab_2d_3d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i  )
566      CALL tab_2d_3d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i  )
567      CALL tab_2d_3d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s  )
568      CALL tab_2d_3d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i )
569      CALL tab_2d_3d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i )
570      CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip )
571      CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip )
572      CALL tab_2d_3d( npti, nptidx(1:npti), v_il_2d(1:npti,1:jpl), v_il )
573      CALL tab_2d_3d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su )
574      DO jl = 1, jpl
575         DO jk = 1, nlay_s
576            CALL tab_1d_2d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) )
577         END DO
578         DO jk = 1, nlay_i
579            CALL tab_1d_2d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
580         END DO
581      END DO
582      !
583   END SUBROUTINE itd_shiftice
584   
585
586   SUBROUTINE ice_itd_reb( kt )
587      !!------------------------------------------------------------------
588      !!                ***  ROUTINE ice_itd_reb ***
589      !!
590      !! ** Purpose : rebin - rebins thicknesses into defined categories
591      !!
592      !! ** Method  : If a category thickness is out of bounds, shift part (for down to top)
593      !!              or entire (for top to down) area, volume, and energy
594      !!              to the neighboring category
595      !!------------------------------------------------------------------
596      INTEGER , INTENT (in) ::   kt      ! Ocean time step
597      INTEGER ::   ji, jj, jl   ! dummy loop indices
598      !
599      INTEGER , DIMENSION(jpij,jpl-1) ::   jdonor           ! donor category index
600      REAL(wp), DIMENSION(jpij,jpl-1) ::   zdaice, zdvice   ! ice area and volume transferred
601      !!------------------------------------------------------------------
602      IF( ln_timing )   CALL timing_start('iceitd_reb')
603      !
604      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_itd_reb: rebining ice thickness distribution' 
605      !
606      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'iceitd_reb', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
607      IF( ln_icediachk )   CALL ice_cons2D  (0, 'iceitd_reb',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft)
608      !
609      jdonor(:,:) = 0
610      zdaice(:,:) = 0._wp
611      zdvice(:,:) = 0._wp
612      !
613      !                       !---------------------------------------
614      DO jl = 1, jpl-1        ! identify thicknesses that are too big
615         !                    !---------------------------------------
616         npti = 0   ;   nptidx(:) = 0
617         DO jj = 1, jpj
618            DO ji = 1, jpi
619               IF( a_i(ji,jj,jl) > 0._wp .AND. v_i(ji,jj,jl) > (a_i(ji,jj,jl) * hi_max(jl)) ) THEN
620                  npti = npti + 1
621                  nptidx( npti ) = (jj - 1) * jpi + ji                 
622               ENDIF
623            END DO
624         END DO
625         !
626         IF( npti > 0 ) THEN           
627            !!clem   CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d(1:npti), h_i(:,:,jl) )
628            CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl) )
629            CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl) )
630            !
631            DO ji = 1, npti
632               jdonor(ji,jl)  = jl 
633               ! how much of a_i you send in cat sup is somewhat arbitrary
634               ! these are from CICE => transfer everything
635               !!zdaice(ji,jl)  = a_i_1d(ji)
636               !!zdvice(ji,jl)  = v_i_1d(ji)
637               ! these are from LLN => transfer only half of the category
638               zdaice(ji,jl)  =                       0.5_wp  * a_i_1d(ji)
639               zdvice(ji,jl)  = v_i_1d(ji) - (1._wp - 0.5_wp) * a_i_1d(ji) * hi_mean(jl)
640            END DO
641            !
642            CALL itd_shiftice( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) )  ! Shift jl=>jl+1
643            ! Reset shift parameters
644            jdonor(1:npti,jl) = 0
645            zdaice(1:npti,jl) = 0._wp
646            zdvice(1:npti,jl) = 0._wp
647         ENDIF
648         !
649      END DO
650
651      !                       !-----------------------------------------
652      DO jl = jpl-1, 1, -1    ! Identify thicknesses that are too small
653         !                    !-----------------------------------------
654         npti = 0 ; nptidx(:) = 0
655         DO jj = 1, jpj
656            DO ji = 1, jpi
657               IF( a_i(ji,jj,jl+1) > 0._wp .AND. v_i(ji,jj,jl+1) <= (a_i(ji,jj,jl+1) * hi_max(jl)) ) THEN
658                  npti = npti + 1
659                  nptidx( npti ) = (jj - 1) * jpi + ji                 
660               ENDIF
661            END DO
662         END DO
663         !
664         IF( npti > 0 ) THEN
665            CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl+1) ) ! jl+1 is ok
666            CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl+1) ) ! jl+1 is ok
667            !
668            DO ji = 1, npti
669               jdonor(ji,jl) = jl + 1
670               zdaice(ji,jl) = a_i_1d(ji) 
671               zdvice(ji,jl) = v_i_1d(ji)
672            END DO
673            !
674            CALL itd_shiftice( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) )  ! Shift jl+1=>jl
675            ! Reset shift parameters
676            jdonor(1:npti,jl) = 0
677            zdaice(1:npti,jl) = 0._wp
678            zdvice(1:npti,jl) = 0._wp
679         ENDIF
680         !
681      END DO
682      !
683      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'iceitd_reb', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
684      IF( ln_icediachk )   CALL ice_cons2D  (1, 'iceitd_reb',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft)
685      IF( ln_timing    )   CALL timing_stop ('iceitd_reb')
686      !
687   END SUBROUTINE ice_itd_reb
688
689
690   SUBROUTINE ice_itd_init
691      !!------------------------------------------------------------------
692      !!                ***  ROUTINE ice_itd_init ***
693      !!
694      !! ** Purpose :   Initializes the ice thickness distribution
695      !! ** Method  :   ...
696      !! ** input   :   Namelist namitd
697      !!-------------------------------------------------------------------
698      INTEGER  ::   jl            ! dummy loop index
699      INTEGER  ::   ios, ioptio   ! Local integer output status for namelist read
700      REAL(wp) ::   zhmax, znum, zden, zalpha   !   -      -
701      !
702      NAMELIST/namitd/ ln_cat_hfn, rn_himean, ln_cat_usr, rn_catbnd, rn_himin, rn_himax
703      !!------------------------------------------------------------------
704      !
705      REWIND( numnam_ice_ref )      ! Namelist namitd in reference namelist : Parameters for ice
706      READ  ( numnam_ice_ref, namitd, IOSTAT = ios, ERR = 901)
707901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namitd in reference namelist' )
708      REWIND( numnam_ice_cfg )      ! Namelist namitd in configuration namelist : Parameters for ice
709      READ  ( numnam_ice_cfg, namitd, IOSTAT = ios, ERR = 902 )
710902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namitd in configuration namelist' )
711      IF(lwm) WRITE( numoni, namitd )
712      !
713      IF(lwp) THEN                  ! control print
714         WRITE(numout,*)
715         WRITE(numout,*) 'ice_itd_init: Initialization of ice cat distribution '
716         WRITE(numout,*) '~~~~~~~~~~~~'
717         WRITE(numout,*) '   Namelist namitd: '
718         WRITE(numout,*) '      Ice categories are defined by a function of rn_himean**(-0.05)    ln_cat_hfn = ', ln_cat_hfn
719         WRITE(numout,*) '         mean ice thickness in the domain                               rn_himean  = ', rn_himean
720         WRITE(numout,*) '      Ice categories are defined by rn_catbnd                           ln_cat_usr = ', ln_cat_usr
721         WRITE(numout,*) '      minimum ice thickness allowed                                     rn_himin   = ', rn_himin 
722         WRITE(numout,*) '      maximum ice thickness allowed                                     rn_himax   = ', rn_himax 
723      ENDIF
724      !
725      !-----------------------------------!
726      !  Thickness categories boundaries  !
727      !-----------------------------------!
728      !                             !== set the choice of ice categories ==!
729      ioptio = 0 
730      IF( ln_cat_hfn ) THEN   ;   ioptio = ioptio + 1   ;   nice_catbnd = np_cathfn    ;   ENDIF
731      IF( ln_cat_usr ) THEN   ;   ioptio = ioptio + 1   ;   nice_catbnd = np_catusr    ;   ENDIF
732      IF( ioptio /= 1 )   CALL ctl_stop( 'ice_itd_init: choose one and only one ice categories boundaries' )
733      !
734      SELECT CASE( nice_catbnd )
735      !                                !------------------------!
736      CASE( np_cathfn )                ! h^(-alpha) function
737         !                             !------------------------!
738         zalpha = 0.05_wp
739         zhmax  = 3._wp * rn_himean
740         hi_max(0) = 0._wp
741         DO jl = 1, jpl
742            znum = jpl * ( zhmax+1 )**zalpha
743            zden = REAL( jpl-jl , wp ) * ( zhmax + 1._wp )**zalpha + REAL( jl , wp )
744            hi_max(jl) = ( znum / zden )**(1./zalpha) - 1
745         END DO
746         !                             !------------------------!
747      CASE( np_catusr )                ! user defined
748         !                             !------------------------!
749         DO jl = 0, jpl
750            hi_max(jl) = rn_catbnd(jl)
751         END DO
752         !
753      END SELECT
754      !
755      DO jl = 1, jpl                ! mean thickness by category
756         hi_mean(jl) = ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp
757      END DO
758      !
759      hi_max(jpl) = rn_himax        ! set to a big value to ensure that all ice is thinner than hi_max(jpl)
760      !
761      IF(lwp) WRITE(numout,*)
762      IF(lwp) WRITE(numout,*) '   ===>>>   resulting thickness category boundaries :'
763      IF(lwp) WRITE(numout,*) '            hi_max(:)= ', hi_max(0:jpl)
764      !
765      IF( hi_max(1) < rn_himin )   CALL ctl_stop('ice_itd_init: the upper bound of the 1st category must be bigger than rn_himin')
766      !
767   END SUBROUTINE ice_itd_init
768
769#else
770   !!----------------------------------------------------------------------
771   !!   Default option :         Empty module         NO SI3 sea-ice model
772   !!----------------------------------------------------------------------
773#endif
774
775   !!======================================================================
776END MODULE iceitd
Note: See TracBrowser for help on using the repository browser.