source: NEMO/trunk/src/ICE/iceitd.F90 @ 13226

Last change on this file since 13226 was 13226, checked in by orioltp, 3 months ago

Merging dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation into the trunk

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