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/trunk/src/ICE – NEMO

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

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

change history of the ice routines

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