New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
iceitd.F90 in NEMO/releases/release-4.0/src/ICE – NEMO

source: NEMO/releases/release-4.0/src/ICE/iceitd.F90 @ 10993

Last change on this file since 10993 was 10993, checked in by clem, 6 years ago

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