source: NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/iceitd.F90 @ 12636

Last change on this file since 12636 was 12636, checked in by dancopsey, 8 months ago

Add:

  • Melt pond lids
  • Melt pond maximum area and thickness
  • Melt pond vertical flushing
  • Area contributing to melt ponds depends on total ice fraction
File size: 37.5 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   !
50   !!----------------------------------------------------------------------
51   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
52   !! $Id$
53   !! Software governed by the CeCILL license (see ./LICENSE)
54   !!----------------------------------------------------------------------
55CONTAINS
56
57   SUBROUTINE ice_itd_rem( kt )
58      !!------------------------------------------------------------------
59      !!                ***  ROUTINE ice_itd_rem ***
60      !!
61      !! ** Purpose :   computes the redistribution of ice thickness
62      !!                after thermodynamic growth of ice thickness
63      !!
64      !! ** Method  :   Linear remapping
65      !!
66      !! References :   W.H. Lipscomb, JGR 2001
67      !!------------------------------------------------------------------
68      INTEGER , INTENT (in) ::   kt      ! Ocean time step
69      !
70      INTEGER  ::   ji, jj, jl, jcat     ! dummy loop index
71      INTEGER  ::   ipti                 ! local integer
72      REAL(wp) ::   zx1, zwk1, zdh0, zetamin, zdamax   ! local scalars
73      REAL(wp) ::   zx2, zwk2, zda0, zetamax           !   -      -
74      REAL(wp) ::   zx3       
75      REAL(wp) ::   zslope          ! used to compute local thermodynamic "speeds"
76      !
77      INTEGER , DIMENSION(jpij)       ::   iptidx          ! compute remapping or not
78      INTEGER , DIMENSION(jpij,jpl-1) ::   jdonor          ! donor category index
79      REAL(wp), DIMENSION(jpij,jpl)   ::   zdhice          ! ice thickness increment
80      REAL(wp), DIMENSION(jpij,jpl)   ::   g0, g1          ! coefficients for fitting the line of the ITD
81      REAL(wp), DIMENSION(jpij,jpl)   ::   hL, hR          ! left and right boundary for the ITD for each thickness
82      REAL(wp), DIMENSION(jpij,jpl-1) ::   zdaice, zdvice  ! local increment of ice area and volume
83      REAL(wp), DIMENSION(jpij)       ::   zhb0, zhb1      ! category boundaries for thinnes categories
84      REAL(wp), DIMENSION(jpij,0:jpl) ::   zhbnew          ! new boundaries of ice categories
85      !!------------------------------------------------------------------
86
87      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_itd_rem: remapping ice thickness distribution' 
88
89      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'iceitd_rem', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
90      IF( ln_icediachk )   CALL ice_cons2D  (0, 'iceitd_rem',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft)
91
92      !-----------------------------------------------------------------------------------------------
93      !  1) Identify grid cells with ice
94      !-----------------------------------------------------------------------------------------------
95      at_i(:,:) = SUM( a_i, dim=3 )
96      !
97      npti = 0   ;   nptidx(:) = 0
98      DO jj = 1, jpj
99         DO ji = 1, jpi
100            IF ( at_i(ji,jj) > epsi10 ) THEN
101               npti = npti + 1
102               nptidx( npti ) = (jj - 1) * jpi + ji
103            ENDIF
104         END DO
105      END DO
106     
107      !-----------------------------------------------------------------------------------------------
108      !  2) Compute new category boundaries
109      !-----------------------------------------------------------------------------------------------
110      IF( npti > 0 ) THEN
111         !
112         zdhice(:,:) = 0._wp
113         zhbnew(:,:) = 0._wp
114         !
115         CALL tab_3d_2d( npti, nptidx(1:npti), h_i_2d (1:npti,1:jpl), h_i   )
116         CALL tab_3d_2d( npti, nptidx(1:npti), h_ib_2d(1:npti,1:jpl), h_i_b )
117         CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i   )
118         CALL tab_3d_2d( npti, nptidx(1:npti), a_ib_2d(1:npti,1:jpl), a_i_b )
119         !
120         DO jl = 1, jpl
121            ! Compute thickness change in each ice category
122            DO ji = 1, npti
123               IF( a_i_2d(ji,jl) > epsi10 )   zdhice(ji,jl) = h_i_2d(ji,jl) - h_ib_2d(ji,jl)
124            END DO
125         END DO
126         !
127         ! --- New boundaries for category 1:jpl-1 --- !
128         DO jl = 1, jpl - 1
129            !
130            DO ji = 1, npti
131               !
132               ! --- New boundary: Hn* = Hn + Fn*dt --- !
133               !     Fn*dt = ( fn + (fn+1 - fn)/(hn+1 - hn) * (Hn - hn) ) * dt = zdhice + zslope * (Hmax - h_i_b)
134               !
135               IF    ( a_ib_2d(ji,jl) >  epsi10 .AND. a_ib_2d(ji,jl+1) >  epsi10 ) THEN   ! a(jl+1) & a(jl) /= 0
136                  zslope        = ( zdhice(ji,jl+1) - zdhice(ji,jl) ) / ( h_ib_2d(ji,jl+1) - h_ib_2d(ji,jl) )
137                  zhbnew(ji,jl) = hi_max(jl) + zdhice(ji,jl) + zslope * ( hi_max(jl) - h_ib_2d(ji,jl) )
138               ELSEIF( a_ib_2d(ji,jl) >  epsi10 .AND. a_ib_2d(ji,jl+1) <= epsi10 ) THEN   ! a(jl+1)=0 => Hn* = Hn + fn*dt
139                  zhbnew(ji,jl) = hi_max(jl) + zdhice(ji,jl)
140               ELSEIF( a_ib_2d(ji,jl) <= epsi10 .AND. a_ib_2d(ji,jl+1) >  epsi10 ) THEN   ! a(jl)=0 => Hn* = Hn + fn+1*dt
141                  zhbnew(ji,jl) = hi_max(jl) + zdhice(ji,jl+1)
142               ELSE                                                                       ! a(jl+1) & a(jl) = 0
143                  zhbnew(ji,jl) = hi_max(jl)
144               ENDIF
145               !
146               ! --- 2 conditions for remapping --- !
147               ! 1) hn(t+1)+espi < Hn* < hn+1(t+1)-epsi               
148               !    Note: hn(t+1) must not be too close to either HR or HL otherwise a division by nearly 0 is possible
149               !          in itd_glinear in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice)
150               IF( a_i_2d(ji,jl  ) > epsi10 .AND. h_i_2d(ji,jl  ) > ( zhbnew(ji,jl) - epsi10 ) )   nptidx(ji) = 0
151               IF( a_i_2d(ji,jl+1) > epsi10 .AND. h_i_2d(ji,jl+1) < ( zhbnew(ji,jl) + epsi10 ) )   nptidx(ji) = 0
152               !
153               ! 2) Hn-1 < Hn* < Hn+1 
154               IF( zhbnew(ji,jl) < hi_max(jl-1) )   nptidx(ji) = 0
155               IF( zhbnew(ji,jl) > hi_max(jl+1) )   nptidx(ji) = 0
156               !
157            END DO
158         END DO
159         !
160         ! --- New boundaries for category jpl --- !
161         DO ji = 1, npti
162            IF( a_i_2d(ji,jpl) > epsi10 ) THEN
163               zhbnew(ji,jpl) = MAX( hi_max(jpl-1), 3._wp * h_i_2d(ji,jpl) - 2._wp * zhbnew(ji,jpl-1) )
164            ELSE
165               zhbnew(ji,jpl) = hi_max(jpl) 
166            ENDIF
167            !
168            ! --- 1 additional condition for remapping (1st category) --- !
169            ! H0+epsi < h1(t) < H1-epsi
170            !    h1(t) must not be too close to either HR or HL otherwise a division by nearly 0 is possible
171            !    in itd_glinear in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice)
172            IF( h_ib_2d(ji,1) < ( hi_max(0) + epsi10 ) )   nptidx(ji) = 0
173            IF( h_ib_2d(ji,1) > ( hi_max(1) - epsi10 ) )   nptidx(ji) = 0
174         END DO
175         !
176         !-----------------------------------------------------------------------------------------------
177         !  3) Identify cells where remapping
178         !-----------------------------------------------------------------------------------------------
179         ipti = 0   ;   iptidx(:) = 0
180         DO ji = 1, npti
181            IF( nptidx(ji) /= 0 ) THEN
182               ipti = ipti + 1
183               iptidx(ipti)   = nptidx(ji)
184               zhbnew(ipti,:) = zhbnew(ji,:)  ! adjust zhbnew to new indices
185            ENDIF
186         END DO
187         nptidx(:) = iptidx(:)
188         npti      = ipti
189         !
190      ENDIF
191   
192      !-----------------------------------------------------------------------------------------------
193      !  4) Compute g(h)
194      !-----------------------------------------------------------------------------------------------
195      IF( npti > 0 ) THEN
196         !
197         zhb0(:) = hi_max(0)   ;   zhb1(:) = hi_max(1)
198         g0(:,:) = 0._wp       ;   g1(:,:) = 0._wp 
199         hL(:,:) = 0._wp       ;   hR(:,:) = 0._wp 
200         !
201         DO jl = 1, jpl
202            !
203            CALL tab_2d_1d( npti, nptidx(1:npti), h_ib_1d(1:npti), h_i_b(:,:,jl) )
204            CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i  (:,:,jl) )
205            CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i  (:,:,jl) )
206            CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d (1:npti), v_i  (:,:,jl) )
207            !
208            IF( jl == 1 ) THEN
209               
210               ! --- g(h) for category 1 --- !
211               CALL itd_glinear( zhb0(1:npti)  , zhb1(1:npti)  , h_ib_1d(1:npti)  , a_i_1d(1:npti)  ,  &   ! in
212                  &              g0  (1:npti,1), g1  (1:npti,1), hL     (1:npti,1), hR    (1:npti,1)   )   ! out
213                  !
214               ! Area lost due to melting of thin ice
215               DO ji = 1, npti
216                  !
217                  IF( a_i_1d(ji) > epsi10 ) THEN
218                     !
219                     zdh0 =  h_i_1d(ji) - h_ib_1d(ji)               
220                     IF( zdh0 < 0.0 ) THEN      !remove area from category 1
221                        zdh0 = MIN( -zdh0, hi_max(1) )
222                        !Integrate g(1) from 0 to dh0 to estimate area melted
223                        zetamax = MIN( zdh0, hR(ji,1) ) - hL(ji,1)
224                        !
225                        IF( zetamax > 0.0 ) THEN
226                           zx1    = zetamax
227                           zx2    = 0.5 * zetamax * zetamax 
228                           zda0   = g1(ji,1) * zx2 + g0(ji,1) * zx1                        ! ice area removed
229                           zdamax = a_i_1d(ji) * (1.0 - h_i_1d(ji) / h_ib_1d(ji) ) ! Constrain new thickness <= h_i               
230                           zda0   = MIN( zda0, zdamax )                                                  ! ice area lost due to melting
231                           !     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_H12 )   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) > 0._wp )  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), lh_ip_2d(1:npti,1:jpl), lh_ip )
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_H12 ) 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                  ztrans          = lh_ip_2d(ji,jl1) * zworka(ji)     ! Pond lid thickness
488                  lh_ip_2d(ji,jl1) = lh_ip_2d(ji,jl1) - ztrans
489                  lh_ip_2d(ji,jl2) = lh_ip_2d(ji,jl2) + ztrans
490               ENDIF
491               !
492            ENDIF   ! jl1 >0
493         END DO
494         !
495         DO jk = 1, nlay_s         !--- Snow heat content
496            DO ji = 1, npti
497               !
498               jl1 = kdonor(ji,jl)
499               !
500               IF( jl1 > 0 ) THEN
501                  IF(jl1 == jl) THEN  ;  jl2 = jl+1
502                  ELSE                ;  jl2 = jl
503                  ENDIF
504                  ztrans             = ze_s_2d(ji,jk,jl1) * zworkv(ji)
505                  ze_s_2d(ji,jk,jl1) = ze_s_2d(ji,jk,jl1) - ztrans
506                  ze_s_2d(ji,jk,jl2) = ze_s_2d(ji,jk,jl2) + ztrans
507               ENDIF
508            END DO
509         END DO
510         !
511         DO jk = 1, nlay_i         !--- Ice heat content
512            DO ji = 1, npti
513               !
514               jl1 = kdonor(ji,jl)
515               !
516               IF( jl1 > 0 ) THEN
517                  IF(jl1 == jl) THEN  ;  jl2 = jl+1
518                  ELSE                ;  jl2 = jl
519                  ENDIF
520                  ztrans             = ze_i_2d(ji,jk,jl1) * zworkv(ji)
521                  ze_i_2d(ji,jk,jl1) = ze_i_2d(ji,jk,jl1) - ztrans
522                  ze_i_2d(ji,jk,jl2) = ze_i_2d(ji,jk,jl2) + ztrans
523               ENDIF
524            END DO
525         END DO
526         !
527      END DO                   ! boundaries, 1 to jpl-1
528
529      !-------------------
530      ! 3) roundoff errors
531      !-------------------
532      ! clem: The transfer between one category to another can lead to very small negative values (-1.e-20)
533      !       because of truncation error ( i.e. 1. - 1. /= 0 )
534      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, lh_ip_2d, ze_s_2d, ze_i_2d )
535
536      ! at_i must be <= rn_amax
537      zworka(1:npti) = SUM( a_i_2d(1:npti,:), dim=2 )
538      DO jl  = 1, jpl
539         WHERE( zworka(1:npti) > rn_amax_1d(1:npti) )   &
540            &   a_i_2d(1:npti,jl) = a_i_2d(1:npti,jl) * rn_amax_1d(1:npti) / zworka(1:npti)
541      END DO
542     
543      !-------------------------------------------------------------------------------
544      ! 4) Update ice thickness and temperature
545      !-------------------------------------------------------------------------------
546      WHERE( a_i_2d(1:npti,:) >= epsi20 )
547         h_i_2d (1:npti,:)  =  v_i_2d(1:npti,:) / a_i_2d(1:npti,:) 
548         t_su_2d(1:npti,:)  =  zaTsfn(1:npti,:) / a_i_2d(1:npti,:) 
549      ELSEWHERE
550         h_i_2d (1:npti,:)  = 0._wp
551         t_su_2d(1:npti,:)  = rt0
552      END WHERE
553      !
554      CALL tab_2d_3d( npti, nptidx(1:npti), h_i_2d (1:npti,1:jpl), h_i  )
555      CALL tab_2d_3d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i  )
556      CALL tab_2d_3d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i  )
557      CALL tab_2d_3d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s  )
558      CALL tab_2d_3d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i )
559      CALL tab_2d_3d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i )
560      CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip )
561      CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip )
562      CALL tab_2d_3d( npti, nptidx(1:npti), lh_ip_2d(1:npti,1:jpl), lh_ip )
563      CALL tab_2d_3d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su )
564      DO jl = 1, jpl
565         DO jk = 1, nlay_s
566            CALL tab_1d_2d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) )
567         END DO
568         DO jk = 1, nlay_i
569            CALL tab_1d_2d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
570         END DO
571      END DO
572      !
573   END SUBROUTINE itd_shiftice
574   
575
576   SUBROUTINE ice_itd_reb( kt )
577      !!------------------------------------------------------------------
578      !!                ***  ROUTINE ice_itd_reb ***
579      !!
580      !! ** Purpose : rebin - rebins thicknesses into defined categories
581      !!
582      !! ** Method  : If a category thickness is out of bounds, shift part (for down to top)
583      !!              or entire (for top to down) area, volume, and energy
584      !!              to the neighboring category
585      !!------------------------------------------------------------------
586      INTEGER , INTENT (in) ::   kt      ! Ocean time step
587      INTEGER ::   ji, jj, jl   ! dummy loop indices
588      !
589      INTEGER , DIMENSION(jpij,jpl-1) ::   jdonor           ! donor category index
590      REAL(wp), DIMENSION(jpij,jpl-1) ::   zdaice, zdvice   ! ice area and volume transferred
591      !!------------------------------------------------------------------
592      !
593      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_itd_reb: rebining ice thickness distribution' 
594      !
595      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'iceitd_reb', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
596      IF( ln_icediachk )   CALL ice_cons2D  (0, 'iceitd_reb',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft)
597      !
598      jdonor(:,:) = 0
599      zdaice(:,:) = 0._wp
600      zdvice(:,:) = 0._wp
601      !
602      !                       !---------------------------------------
603      DO jl = 1, jpl-1        ! identify thicknesses that are too big
604         !                    !---------------------------------------
605         npti = 0   ;   nptidx(:) = 0
606         DO jj = 1, jpj
607            DO ji = 1, jpi
608               IF( a_i(ji,jj,jl) > 0._wp .AND. v_i(ji,jj,jl) > (a_i(ji,jj,jl) * hi_max(jl)) ) THEN
609                  npti = npti + 1
610                  nptidx( npti ) = (jj - 1) * jpi + ji                 
611               ENDIF
612            END DO
613         END DO
614         !
615!!clem   CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d(1:npti), h_i(:,:,jl) )
616         CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl) )
617         CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl) )
618         !
619         DO ji = 1, npti
620            jdonor(ji,jl)  = jl 
621            ! how much of a_i you send in cat sup is somewhat arbitrary
622!!clem: these do not work properly after a restart (I do not know why) => not sure it is still true
623!!          zdaice(ji,jl)  = a_i_1d(ji) * ( h_i_1d(ji) - hi_max(jl) + epsi10 ) / h_i_1d(ji) 
624!!          zdvice(ji,jl)  = v_i_1d(ji) - ( a_i_1d(ji) - zdaice(ji,jl) ) * ( hi_max(jl) - epsi10 )
625!!clem: these do not work properly after a restart (I do not know why) => not sure it is still true
626!!          zdaice(ji,jl)  = a_i_1d(ji)
627!!          zdvice(ji,jl)  = v_i_1d(ji)
628!!clem: these are from UCL and work ok
629            zdaice(ji,jl)  = a_i_1d(ji) * 0.5_wp
630            zdvice(ji,jl)  = v_i_1d(ji) - zdaice(ji,jl) * ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp
631         END DO
632         !
633         IF( npti > 0 ) THEN
634            CALL itd_shiftice( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) )  ! Shift jl=>jl+1
635            ! Reset shift parameters
636            jdonor(1:npti,jl) = 0
637            zdaice(1:npti,jl) = 0._wp
638            zdvice(1:npti,jl) = 0._wp
639         ENDIF
640         !
641      END DO
642
643      !                       !-----------------------------------------
644      DO jl = jpl-1, 1, -1    ! Identify thicknesses that are too small
645         !                    !-----------------------------------------
646         npti = 0 ; nptidx(:) = 0
647         DO jj = 1, jpj
648            DO ji = 1, jpi
649               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
650                  npti = npti + 1
651                  nptidx( npti ) = (jj - 1) * jpi + ji                 
652               ENDIF
653            END DO
654         END DO
655         !
656         CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl+1) ) ! jl+1 is ok
657         CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl+1) ) ! jl+1 is ok
658         DO ji = 1, npti
659            jdonor(ji,jl) = jl + 1
660            zdaice(ji,jl) = a_i_1d(ji) 
661            zdvice(ji,jl) = v_i_1d(ji)
662         END DO
663         !
664         IF( npti > 0 ) THEN
665            CALL itd_shiftice( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) )  ! Shift jl+1=>jl
666            ! Reset shift parameters
667            jdonor(1:npti,jl) = 0
668            zdaice(1:npti,jl) = 0._wp
669            zdvice(1:npti,jl) = 0._wp
670         ENDIF
671         !
672      END DO
673      !
674      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'iceitd_reb', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
675      IF( ln_icediachk )   CALL ice_cons2D  (1, 'iceitd_reb',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft)
676      !
677   END SUBROUTINE ice_itd_reb
678
679
680   SUBROUTINE ice_itd_init
681      !!------------------------------------------------------------------
682      !!                ***  ROUTINE ice_itd_init ***
683      !!
684      !! ** Purpose :   Initializes the ice thickness distribution
685      !! ** Method  :   ...
686      !! ** input   :   Namelist namitd
687      !!-------------------------------------------------------------------
688      INTEGER  ::   jl            ! dummy loop index
689      INTEGER  ::   ios, ioptio   ! Local integer output status for namelist read
690      REAL(wp) ::   zhmax, znum, zden, zalpha   !   -      -
691      !
692      NAMELIST/namitd/ ln_cat_hfn, rn_himean, ln_cat_usr, rn_catbnd, rn_himin
693      !!------------------------------------------------------------------
694      !
695      REWIND( numnam_ice_ref )      ! Namelist namitd in reference namelist : Parameters for ice
696      READ  ( numnam_ice_ref, namitd, IOSTAT = ios, ERR = 901)
697901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namitd in reference namelist' )
698      REWIND( numnam_ice_cfg )      ! Namelist namitd in configuration namelist : Parameters for ice
699      READ  ( numnam_ice_cfg, namitd, IOSTAT = ios, ERR = 902 )
700902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namitd in configuration namelist' )
701      IF(lwm) WRITE( numoni, namitd )
702      !
703      IF(lwp) THEN                  ! control print
704         WRITE(numout,*)
705         WRITE(numout,*) 'ice_itd_init: Initialization of ice cat distribution '
706         WRITE(numout,*) '~~~~~~~~~~~~'
707         WRITE(numout,*) '   Namelist namitd: '
708         WRITE(numout,*) '      Ice categories are defined by a function of rn_himean**(-0.05)    ln_cat_hfn = ', ln_cat_hfn
709         WRITE(numout,*) '         mean ice thickness in the domain                               rn_himean  = ', rn_himean
710         WRITE(numout,*) '      Ice categories are defined by rn_catbnd                           ln_cat_usr = ', ln_cat_usr
711         WRITE(numout,*) '      minimum ice thickness                                             rn_himin   = ', rn_himin 
712      ENDIF
713      !
714      !-----------------------------------!
715      !  Thickness categories boundaries  !
716      !-----------------------------------!
717      !                             !== set the choice of ice categories ==!
718      ioptio = 0 
719      IF( ln_cat_hfn ) THEN   ;   ioptio = ioptio + 1   ;   nice_catbnd = np_cathfn    ;   ENDIF
720      IF( ln_cat_usr ) THEN   ;   ioptio = ioptio + 1   ;   nice_catbnd = np_catusr    ;   ENDIF
721      IF( ioptio /= 1 )   CALL ctl_stop( 'ice_itd_init: choose one and only one ice categories boundaries' )
722      !
723      SELECT CASE( nice_catbnd )
724      !                                !------------------------!
725      CASE( np_cathfn )                ! h^(-alpha) function
726         !                             !------------------------!
727         zalpha = 0.05_wp
728         zhmax  = 3._wp * rn_himean
729         hi_max(0) = 0._wp
730         DO jl = 1, jpl
731            znum = jpl * ( zhmax+1 )**zalpha
732            zden = REAL( jpl-jl , wp ) * ( zhmax + 1._wp )**zalpha + REAL( jl , wp )
733            hi_max(jl) = ( znum / zden )**(1./zalpha) - 1
734         END DO
735         !                             !------------------------!
736      CASE( np_catusr )                ! user defined
737         !                             !------------------------!
738         DO jl = 0, jpl
739            hi_max(jl) = rn_catbnd(jl)
740         END DO
741         !
742      END SELECT
743      !
744      DO jl = 1, jpl                ! mean thickness by category
745         hi_mean(jl) = ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp
746      END DO
747      !
748      hi_max(jpl) = 99._wp          ! set to a big value to ensure that all ice is thinner than hi_max(jpl)
749      !
750      IF(lwp) WRITE(numout,*)
751      IF(lwp) WRITE(numout,*) '   ===>>>   resulting thickness category boundaries :'
752      IF(lwp) WRITE(numout,*) '            hi_max(:)= ', hi_max(0:jpl)
753      !
754      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')
755      !
756   END SUBROUTINE ice_itd_init
757
758#else
759   !!----------------------------------------------------------------------
760   !!   Default option :         Empty module         NO SI3 sea-ice model
761   !!----------------------------------------------------------------------
762#endif
763
764   !!======================================================================
765END MODULE iceitd
Note: See TracBrowser for help on using the repository browser.