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

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

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

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

r4_trunk: first change of DO loops for routines to be merged, see #2523
  • Property svn:keywords set to Id
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   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_2D_11_11
100         IF ( at_i(ji,jj) > epsi10 ) THEN
101            npti = npti + 1
102            nptidx( npti ) = (jj - 1) * jpi + ji
103         ENDIF
104      END_2D
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 of thin ice (zdamax > 0)
230                           ! Remove area, conserving volume
231                           h_i_1d(ji) = h_i_1d(ji) * a_i_1d(ji) / ( a_i_1d(ji) - zda0 )
232                           a_i_1d(ji) = a_i_1d(ji) - zda0
233                           v_i_1d(ji) = a_i_1d(ji) * h_i_1d(ji) ! useless ?
234                        ENDIF
235                        !
236                     ELSE ! if ice accretion zdh0 > 0
237                        ! zhbnew was 0, and is shifted to the right to account for thin ice growth in openwater (F0 = f1)
238                        zhbnew(ji,0) = MIN( zdh0, hi_max(1) ) 
239                     ENDIF
240                     !
241                  ENDIF
242                  !
243               END DO
244               !
245               CALL tab_1d_2d( npti, nptidx(1:npti), h_i_1d(1:npti), h_i(:,:,jl) )
246               CALL tab_1d_2d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl) )
247               CALL tab_1d_2d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl) )
248               !
249            ENDIF ! jl=1
250            !
251            ! --- g(h) for each thickness category --- ! 
252            CALL itd_glinear( zhbnew(1:npti,jl-1), zhbnew(1:npti,jl), h_i_1d(1:npti)   , a_i_1d(1:npti)   ,  &   ! in
253               &              g0    (1:npti,jl  ), g1    (1:npti,jl), hL    (1:npti,jl), hR    (1:npti,jl)   )   ! out
254            !
255         END DO
256         
257         !-----------------------------------------------------------------------------------------------
258         !  5) Compute area and volume to be shifted across each boundary (Eq. 18)
259         !-----------------------------------------------------------------------------------------------
260         DO jl = 1, jpl - 1
261            !
262            DO ji = 1, npti
263               !
264               ! left and right integration limits in eta space
265               IF (zhbnew(ji,jl) > hi_max(jl)) THEN ! Hn* > Hn => transfer from jl to jl+1
266                  zetamin = MAX( hi_max(jl)   , hL(ji,jl) ) - hL(ji,jl)   ! hi_max(jl) - hL
267                  zetamax = MIN( zhbnew(ji,jl), hR(ji,jl) ) - hL(ji,jl)   ! hR - hL
268                  jdonor(ji,jl) = jl
269               ELSE                                 ! Hn* <= Hn => transfer from jl+1 to jl
270                  zetamin = 0.0
271                  zetamax = MIN( hi_max(jl), hR(ji,jl+1) ) - hL(ji,jl+1)  ! hi_max(jl) - hL
272                  jdonor(ji,jl) = jl + 1
273               ENDIF
274               zetamax = MAX( zetamax, zetamin ) ! no transfer if etamax < etamin
275               !
276               zx1  = zetamax - zetamin
277               zwk1 = zetamin * zetamin
278               zwk2 = zetamax * zetamax
279               zx2  = 0.5 * ( zwk2 - zwk1 )
280               zwk1 = zwk1 * zetamin
281               zwk2 = zwk2 * zetamax
282               zx3  = 1.0 / 3.0 * ( zwk2 - zwk1 )
283               jcat = jdonor(ji,jl)
284               zdaice(ji,jl) = g1(ji,jcat)*zx2 + g0(ji,jcat)*zx1
285               zdvice(ji,jl) = g1(ji,jcat)*zx3 + g0(ji,jcat)*zx2 + zdaice(ji,jl)*hL(ji,jcat)
286               !
287            END DO
288         END DO
289         
290         !----------------------------------------------------------------------------------------------
291         ! 6) Shift ice between categories
292         !----------------------------------------------------------------------------------------------
293         CALL itd_shiftice ( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) )
294         
295         !----------------------------------------------------------------------------------------------
296         ! 7) Make sure h_i >= minimum ice thickness hi_min
297         !----------------------------------------------------------------------------------------------
298         CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i (:,:,1) )
299         CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i (:,:,1) )
300         CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d(1:npti), a_ip(:,:,1) )
301         !
302         DO ji = 1, npti
303            IF ( a_i_1d(ji) > epsi10 .AND. h_i_1d(ji) < rn_himin ) THEN
304               a_i_1d(ji) = a_i_1d(ji) * h_i_1d(ji) / rn_himin 
305               IF( ln_pnd_LEV )   a_ip_1d(ji) = a_ip_1d(ji) * h_i_1d(ji) / rn_himin
306               h_i_1d(ji) = rn_himin
307            ENDIF
308         END DO
309         !
310         CALL tab_1d_2d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i (:,:,1) )
311         CALL tab_1d_2d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i (:,:,1) )
312         CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d(1:npti), a_ip(:,:,1) )
313         !
314      ENDIF
315      !
316      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'iceitd_rem', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
317      IF( ln_icediachk )   CALL ice_cons2D  (1, 'iceitd_rem',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_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) > epsi10 )  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), v_il_2d(1:npti,1:jpl), v_il )
412      CALL tab_3d_2d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su )
413      DO jl = 1, jpl
414         DO jk = 1, nlay_s
415            CALL tab_2d_1d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) )
416         END DO
417         DO jk = 1, nlay_i
418            CALL tab_2d_1d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
419         END DO
420      END DO
421      ! to correct roundoff errors on a_i
422      CALL tab_2d_1d( npti, nptidx(1:npti), rn_amax_1d(1:npti), rn_amax_2d )
423
424      !----------------------------------------------------------------------------------------------
425      ! 1) Define a variable equal to a_i*T_su
426      !----------------------------------------------------------------------------------------------
427      DO jl = 1, jpl
428         DO ji = 1, npti
429            zaTsfn(ji,jl) = a_i_2d(ji,jl) * t_su_2d(ji,jl)
430         END DO
431      END DO
432     
433      !-------------------------------------------------------------------------------
434      ! 2) Transfer volume and energy between categories
435      !-------------------------------------------------------------------------------
436      DO jl = 1, jpl - 1
437         DO ji = 1, npti
438            !
439            jl1 = kdonor(ji,jl)
440            !
441            IF( jl1 > 0 ) THEN
442               !
443               IF ( jl1 == jl  ) THEN   ;   jl2 = jl1+1
444               ELSE                     ;   jl2 = jl 
445               ENDIF
446               !
447               IF( v_i_2d(ji,jl1) >= epsi10 ) THEN   ;   zworkv(ji) = pdvice(ji,jl) / v_i_2d(ji,jl1)
448               ELSE                                  ;   zworkv(ji) = 0._wp
449               ENDIF
450               IF( a_i_2d(ji,jl1) >= epsi10 ) THEN   ;   zworka(ji) = pdaice(ji,jl) / a_i_2d(ji,jl1)
451               ELSE                                  ;   zworka(ji) = 0._wp
452               ENDIF
453               !
454               a_i_2d(ji,jl1) = a_i_2d(ji,jl1) - pdaice(ji,jl)       ! Ice areas
455               a_i_2d(ji,jl2) = a_i_2d(ji,jl2) + pdaice(ji,jl)
456               !
457               v_i_2d(ji,jl1) = v_i_2d(ji,jl1) - pdvice(ji,jl)       ! Ice volumes
458               v_i_2d(ji,jl2) = v_i_2d(ji,jl2) + pdvice(ji,jl)
459               !
460               ztrans         = v_s_2d(ji,jl1) * zworkv(ji)          ! Snow volumes
461               v_s_2d(ji,jl1) = v_s_2d(ji,jl1) - ztrans
462               v_s_2d(ji,jl2) = v_s_2d(ji,jl2) + ztrans 
463               !
464               ztrans          = oa_i_2d(ji,jl1) * zworka(ji)        ! Ice age
465               oa_i_2d(ji,jl1) = oa_i_2d(ji,jl1) - ztrans
466               oa_i_2d(ji,jl2) = oa_i_2d(ji,jl2) + ztrans
467               !
468               ztrans          = sv_i_2d(ji,jl1) * zworkv(ji)        ! Ice salinity
469               sv_i_2d(ji,jl1) = sv_i_2d(ji,jl1) - ztrans
470               sv_i_2d(ji,jl2) = sv_i_2d(ji,jl2) + ztrans
471               !
472               ztrans          = zaTsfn(ji,jl1) * zworka(ji)         ! Surface temperature
473               zaTsfn(ji,jl1)  = zaTsfn(ji,jl1) - ztrans
474               zaTsfn(ji,jl2)  = zaTsfn(ji,jl2) + ztrans
475               ! 
476               IF ( ln_pnd_LEV ) THEN
477                  ztrans          = a_ip_2d(ji,jl1) * zworka(ji)     ! Pond fraction
478                  a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - ztrans
479                  a_ip_2d(ji,jl2) = a_ip_2d(ji,jl2) + ztrans
480                  !                                             
481                  ztrans          = v_ip_2d(ji,jl1) * zworka(ji)     ! Pond volume (also proportional to da/a)
482                  v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - ztrans
483                  v_ip_2d(ji,jl2) = v_ip_2d(ji,jl2) + ztrans
484                  !
485                  IF ( ln_pnd_lids ) THEN                            ! Pond lid volume
486                     ztrans          = v_il_2d(ji,jl1) * zworka(ji)
487                     v_il_2d(ji,jl1) = v_il_2d(ji,jl1) - ztrans
488                     v_il_2d(ji,jl2) = v_il_2d(ji,jl2) + ztrans
489                  ENDIF
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, v_il_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), v_il_2d(1:npti,1:jpl), v_il )
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_2D_11_11
607            IF( a_i(ji,jj,jl) > 0._wp .AND. v_i(ji,jj,jl) > (a_i(ji,jj,jl) * hi_max(jl)) ) THEN
608               npti = npti + 1
609               nptidx( npti ) = (jj - 1) * jpi + ji                 
610            ENDIF
611         END_2D
612         !
613!!clem   CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d(1:npti), h_i(:,:,jl) )
614         CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl) )
615         CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl) )
616         !
617         DO ji = 1, npti
618            jdonor(ji,jl)  = jl 
619            ! how much of a_i you send in cat sup is somewhat arbitrary
620!!clem: these do not work properly after a restart (I do not know why) => not sure it is still true
621!!          zdaice(ji,jl)  = a_i_1d(ji) * ( h_i_1d(ji) - hi_max(jl) + epsi10 ) / h_i_1d(ji) 
622!!          zdvice(ji,jl)  = v_i_1d(ji) - ( a_i_1d(ji) - zdaice(ji,jl) ) * ( hi_max(jl) - epsi10 )
623!!clem: these do not work properly after a restart (I do not know why) => not sure it is still true
624!!          zdaice(ji,jl)  = a_i_1d(ji)
625!!          zdvice(ji,jl)  = v_i_1d(ji)
626!!clem: these are from UCL and work ok
627            zdaice(ji,jl)  = a_i_1d(ji) * 0.5_wp
628            zdvice(ji,jl)  = v_i_1d(ji) - zdaice(ji,jl) * ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp
629         END DO
630         !
631         IF( npti > 0 ) THEN
632            CALL itd_shiftice( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) )  ! Shift jl=>jl+1
633            ! Reset shift parameters
634            jdonor(1:npti,jl) = 0
635            zdaice(1:npti,jl) = 0._wp
636            zdvice(1:npti,jl) = 0._wp
637         ENDIF
638         !
639      END DO
640
641      !                       !-----------------------------------------
642      DO jl = jpl-1, 1, -1    ! Identify thicknesses that are too small
643         !                    !-----------------------------------------
644         npti = 0 ; nptidx(:) = 0
645         DO_2D_11_11
646            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
647               npti = npti + 1
648               nptidx( npti ) = (jj - 1) * jpi + ji                 
649            ENDIF
650         END_2D
651         !
652         CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl+1) ) ! jl+1 is ok
653         CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl+1) ) ! jl+1 is ok
654         DO ji = 1, npti
655            jdonor(ji,jl) = jl + 1
656            zdaice(ji,jl) = a_i_1d(ji) 
657            zdvice(ji,jl) = v_i_1d(ji)
658         END DO
659         !
660         IF( npti > 0 ) THEN
661            CALL itd_shiftice( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) )  ! Shift jl+1=>jl
662            ! Reset shift parameters
663            jdonor(1:npti,jl) = 0
664            zdaice(1:npti,jl) = 0._wp
665            zdvice(1:npti,jl) = 0._wp
666         ENDIF
667         !
668      END DO
669      !
670      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'iceitd_reb', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
671      IF( ln_icediachk )   CALL ice_cons2D  (1, 'iceitd_reb',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft)
672      !
673   END SUBROUTINE ice_itd_reb
674
675
676   SUBROUTINE ice_itd_init
677      !!------------------------------------------------------------------
678      !!                ***  ROUTINE ice_itd_init ***
679      !!
680      !! ** Purpose :   Initializes the ice thickness distribution
681      !! ** Method  :   ...
682      !! ** input   :   Namelist namitd
683      !!-------------------------------------------------------------------
684      INTEGER  ::   jl            ! dummy loop index
685      INTEGER  ::   ios, ioptio   ! Local integer output status for namelist read
686      REAL(wp) ::   zhmax, znum, zden, zalpha   !   -      -
687      !
688      NAMELIST/namitd/ ln_cat_hfn, rn_himean, ln_cat_usr, rn_catbnd, rn_himin, rn_himax
689      !!------------------------------------------------------------------
690      !
691      REWIND( numnam_ice_ref )      ! Namelist namitd in reference namelist : Parameters for ice
692      READ  ( numnam_ice_ref, namitd, IOSTAT = ios, ERR = 901)
693901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namitd in reference namelist' )
694      REWIND( numnam_ice_cfg )      ! Namelist namitd in configuration namelist : Parameters for ice
695      READ  ( numnam_ice_cfg, namitd, IOSTAT = ios, ERR = 902 )
696902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namitd in configuration namelist' )
697      IF(lwm) WRITE( numoni, namitd )
698      !
699      IF(lwp) THEN                  ! control print
700         WRITE(numout,*)
701         WRITE(numout,*) 'ice_itd_init: Initialization of ice cat distribution '
702         WRITE(numout,*) '~~~~~~~~~~~~'
703         WRITE(numout,*) '   Namelist namitd: '
704         WRITE(numout,*) '      Ice categories are defined by a function of rn_himean**(-0.05)    ln_cat_hfn = ', ln_cat_hfn
705         WRITE(numout,*) '         mean ice thickness in the domain                               rn_himean  = ', rn_himean
706         WRITE(numout,*) '      Ice categories are defined by rn_catbnd                           ln_cat_usr = ', ln_cat_usr
707         WRITE(numout,*) '      minimum ice thickness allowed                                     rn_himin   = ', rn_himin 
708         WRITE(numout,*) '      maximum ice thickness allowed                                     rn_himax   = ', rn_himax 
709      ENDIF
710      !
711      !-----------------------------------!
712      !  Thickness categories boundaries  !
713      !-----------------------------------!
714      !                             !== set the choice of ice categories ==!
715      ioptio = 0 
716      IF( ln_cat_hfn ) THEN   ;   ioptio = ioptio + 1   ;   nice_catbnd = np_cathfn    ;   ENDIF
717      IF( ln_cat_usr ) THEN   ;   ioptio = ioptio + 1   ;   nice_catbnd = np_catusr    ;   ENDIF
718      IF( ioptio /= 1 )   CALL ctl_stop( 'ice_itd_init: choose one and only one ice categories boundaries' )
719      !
720      SELECT CASE( nice_catbnd )
721      !                                !------------------------!
722      CASE( np_cathfn )                ! h^(-alpha) function
723         !                             !------------------------!
724         zalpha = 0.05_wp
725         zhmax  = 3._wp * rn_himean
726         hi_max(0) = 0._wp
727         DO jl = 1, jpl
728            znum = jpl * ( zhmax+1 )**zalpha
729            zden = REAL( jpl-jl , wp ) * ( zhmax + 1._wp )**zalpha + REAL( jl , wp )
730            hi_max(jl) = ( znum / zden )**(1./zalpha) - 1
731         END DO
732         !                             !------------------------!
733      CASE( np_catusr )                ! user defined
734         !                             !------------------------!
735         DO jl = 0, jpl
736            hi_max(jl) = rn_catbnd(jl)
737         END DO
738         !
739      END SELECT
740      !
741      DO jl = 1, jpl                ! mean thickness by category
742         hi_mean(jl) = ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp
743      END DO
744      !
745      hi_max(jpl) = rn_himax        ! set to a big value to ensure that all ice is thinner than hi_max(jpl)
746      !
747      IF(lwp) WRITE(numout,*)
748      IF(lwp) WRITE(numout,*) '   ===>>>   resulting thickness category boundaries :'
749      IF(lwp) WRITE(numout,*) '            hi_max(:)= ', hi_max(0:jpl)
750      !
751      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')
752      !
753   END SUBROUTINE ice_itd_init
754
755#else
756   !!----------------------------------------------------------------------
757   !!   Default option :         Empty module         NO SI3 sea-ice model
758   !!----------------------------------------------------------------------
759#endif
760
761   !!======================================================================
762END MODULE iceitd
Note: See TracBrowser for help on using the repository browser.