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/branches/2020/dev_r2052_ENHANCE-09_rbourdal_massfluxconvection/src/ICE – NEMO

source: NEMO/branches/2020/dev_r2052_ENHANCE-09_rbourdal_massfluxconvection/src/ICE/iceitd.F90 @ 14009

Last change on this file since 14009 was 14009, checked in by gsamson, 3 years ago

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