source: NEMO/branches/2020/temporary_r4_trunk/src/ICE/iceitd.F90 @ 13466

Last change on this file since 13466 was 13466, checked in by smasson, 7 months ago

r4_trunk: merge r4 13280:13310, see #2523

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