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/trunk/src/ICE – NEMO

source: NEMO/trunk/src/ICE/iceitd.F90 @ 13472

Last change on this file since 13472 was 13472, checked in by smasson, 4 years ago

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