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

source: NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/ICE/iceitd.F90 @ 9957

Last change on this file since 9957 was 9957, checked in by acc, 6 years ago

New development branch for the adaptive-implicit vertical advection (05_Gurvan-Vertical_advection)

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)
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               ! clem: The transfer between one category to another can lead to very small negative values (-1.e-20)
439               !       because of truncation error ( i.e. 1. - 1. /= 0 )
440               !       I do not think it should be a concern since small areas and volumes are erased (in ice_var_zapsmall.F90)
441               !
442               a_i_2d(ji,jl1) = a_i_2d(ji,jl1) - pdaice(ji,jl)       ! Ice areas
443               a_i_2d(ji,jl2) = a_i_2d(ji,jl2) + pdaice(ji,jl)
444               !
445               v_i_2d(ji,jl1) = v_i_2d(ji,jl1) - pdvice(ji,jl)       ! Ice volumes
446               v_i_2d(ji,jl2) = v_i_2d(ji,jl2) + pdvice(ji,jl)
447               !
448               ztrans         = v_s_2d(ji,jl1) * zworkv(ji)          ! Snow volumes
449               v_s_2d(ji,jl1) = v_s_2d(ji,jl1) - ztrans
450               v_s_2d(ji,jl2) = v_s_2d(ji,jl2) + ztrans 
451               !
452               ztrans          = oa_i_2d(ji,jl1) * zworka(ji)        ! Ice age
453               oa_i_2d(ji,jl1) = oa_i_2d(ji,jl1) - ztrans
454               oa_i_2d(ji,jl2) = oa_i_2d(ji,jl2) + ztrans
455               !
456               ztrans          = sv_i_2d(ji,jl1) * zworkv(ji)        ! Ice salinity
457               sv_i_2d(ji,jl1) = sv_i_2d(ji,jl1) - ztrans
458               sv_i_2d(ji,jl2) = sv_i_2d(ji,jl2) + ztrans
459               !
460               ztrans          = zaTsfn(ji,jl1) * zworka(ji)         ! Surface temperature
461               zaTsfn(ji,jl1)  = zaTsfn(ji,jl1) - ztrans
462               zaTsfn(ji,jl2)  = zaTsfn(ji,jl2) + ztrans
463               
464               IF ( ln_pnd_H12 ) THEN
465                  ztrans          = a_ip_2d(ji,jl1) * zworka(ji)     ! Pond fraction
466                  a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - ztrans
467                  a_ip_2d(ji,jl2) = a_ip_2d(ji,jl2) + ztrans
468                  !                                             
469                  ztrans          = v_ip_2d(ji,jl1) * zworka(ji)     ! Pond volume (also proportional to da/a)
470                  v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - ztrans
471                  v_ip_2d(ji,jl2) = v_ip_2d(ji,jl2) + ztrans
472               ENDIF
473               !
474            ENDIF   ! jl1 >0
475         END DO
476         !
477         DO jk = 1, nlay_s         !--- Snow heat content
478            !
479            DO ji = 1, npti
480               ii = MOD( nptidx(ji) - 1, jpi ) + 1
481               ij = ( nptidx(ji) - 1 ) / jpi + 1
482               !
483               jl1 = kdonor(ji,jl)
484               !
485               IF( jl1 > 0 ) THEN
486                  IF(jl1 == jl) THEN  ;  jl2 = jl+1
487                  ELSE                ;  jl2 = jl
488                  ENDIF
489                  !
490                  ztrans            = e_s(ii,ij,jk,jl1) * zworkv(ji)
491                  e_s(ii,ij,jk,jl1) = e_s(ii,ij,jk,jl1) - ztrans
492                  e_s(ii,ij,jk,jl2) = e_s(ii,ij,jk,jl2) + ztrans
493               ENDIF
494            END DO
495         END DO
496         !
497         DO jk = 1, nlay_i         !--- Ice heat content
498            DO ji = 1, npti
499               ii = MOD( nptidx(ji) - 1, jpi ) + 1
500               ij = ( nptidx(ji) - 1 ) / jpi + 1
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                  !
509                  ztrans            = e_i(ii,ij,jk,jl1) * zworkv(ji)
510                  e_i(ii,ij,jk,jl1) = e_i(ii,ij,jk,jl1) - ztrans
511                  e_i(ii,ij,jk,jl2) = e_i(ii,ij,jk,jl2) + ztrans
512               ENDIF
513            END DO
514         END DO
515         !
516      END DO                   ! boundaries, 1 to jpl-1
517     
518      !-------------------------------------------------------------------------------
519      ! 3) Update ice thickness and temperature
520      !-------------------------------------------------------------------------------
521      WHERE( a_i_2d(1:npti,:) >= epsi20 )
522         h_i_2d (1:npti,:)  =  v_i_2d(1:npti,:) / a_i_2d(1:npti,:) 
523         t_su_2d(1:npti,:)  =  zaTsfn(1:npti,:) / a_i_2d(1:npti,:) 
524      ELSEWHERE
525         h_i_2d (1:npti,:)  = 0._wp
526         t_su_2d(1:npti,:)  = rt0
527      END WHERE
528      !
529      CALL tab_2d_3d( npti, nptidx(1:npti), h_i_2d (1:npti,1:jpl), h_i  )
530      CALL tab_2d_3d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i  )
531      CALL tab_2d_3d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i  )
532      CALL tab_2d_3d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s  )
533      CALL tab_2d_3d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i )
534      CALL tab_2d_3d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i )
535      CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip )
536      CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip )
537      CALL tab_2d_3d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su )
538      !
539   END SUBROUTINE itd_shiftice
540   
541
542   SUBROUTINE ice_itd_reb( kt )
543      !!------------------------------------------------------------------
544      !!                ***  ROUTINE ice_itd_reb ***
545      !!
546      !! ** Purpose : rebin - rebins thicknesses into defined categories
547      !!
548      !! ** Method  : If a category thickness is out of bounds, shift part (for down to top)
549      !!              or entire (for top to down) area, volume, and energy
550      !!              to the neighboring category
551      !!------------------------------------------------------------------
552      INTEGER , INTENT (in) ::   kt      ! Ocean time step
553      INTEGER ::   ji, jj, jl   ! dummy loop indices
554      !
555      INTEGER , DIMENSION(jpij,jpl-1) ::   jdonor           ! donor category index
556      REAL(wp), DIMENSION(jpij,jpl-1) ::   zdaice, zdvice   ! ice area and volume transferred
557      !!------------------------------------------------------------------
558      !
559      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_itd_reb: rebining ice thickness distribution' 
560      !
561      jdonor(:,:) = 0
562      zdaice(:,:) = 0._wp
563      zdvice(:,:) = 0._wp
564      !
565      !                       !---------------------------------------
566      DO jl = 1, jpl-1        ! identify thicknesses that are too big
567         !                    !---------------------------------------
568         npti = 0   ;   nptidx(:) = 0
569         DO jj = 1, jpj
570            DO ji = 1, jpi
571               IF( a_i(ji,jj,jl) > 0._wp .AND. v_i(ji,jj,jl) > (a_i(ji,jj,jl) * hi_max(jl)) ) THEN
572                  npti = npti + 1
573                  nptidx( npti ) = (jj - 1) * jpi + ji                 
574               ENDIF
575            END DO
576         END DO
577         !
578!!clem   CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d(1:npti), h_i(:,:,jl) )
579         CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl) )
580         CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl) )
581         !
582         DO ji = 1, npti
583            jdonor(ji,jl)  = jl 
584            ! how much of a_i you send in cat sup is somewhat arbitrary
585!!clem: these do not work properly after a restart (I do not know why) => not sure it is still true
586!!          zdaice(ji,jl)  = a_i_1d(ji) * ( h_i_1d(ji) - hi_max(jl) + epsi10 ) / h_i_1d(ji) 
587!!          zdvice(ji,jl)  = v_i_1d(ji) - ( a_i_1d(ji) - zdaice(ji,jl) ) * ( hi_max(jl) - epsi10 )
588!!clem: these do not work properly after a restart (I do not know why) => not sure it is still true
589!!          zdaice(ji,jl)  = a_i_1d(ji)
590!!          zdvice(ji,jl)  = v_i_1d(ji)
591!!clem: these are from UCL and work ok
592            zdaice(ji,jl)  = a_i_1d(ji) * 0.5_wp
593            zdvice(ji,jl)  = v_i_1d(ji) - zdaice(ji,jl) * ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp
594         END DO
595         !
596         IF( npti > 0 ) THEN
597            CALL itd_shiftice( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) )  ! Shift jl=>jl+1
598            ! Reset shift parameters
599            jdonor(1:npti,jl) = 0
600            zdaice(1:npti,jl) = 0._wp
601            zdvice(1:npti,jl) = 0._wp
602         ENDIF
603         !
604      END DO
605
606      !                       !-----------------------------------------
607      DO jl = jpl-1, 1, -1    ! Identify thicknesses that are too small
608         !                    !-----------------------------------------
609         npti = 0 ; nptidx(:) = 0
610         DO jj = 1, jpj
611            DO ji = 1, jpi
612               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
613                  npti = npti + 1
614                  nptidx( npti ) = (jj - 1) * jpi + ji                 
615               ENDIF
616            END DO
617         END DO
618         !
619         CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl+1) ) ! jl+1 is ok
620         CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl+1) ) ! jl+1 is ok
621         DO ji = 1, npti
622            jdonor(ji,jl) = jl + 1
623            zdaice(ji,jl) = a_i_1d(ji) 
624            zdvice(ji,jl) = v_i_1d(ji)
625         END DO
626         !
627         IF( npti > 0 ) THEN
628            CALL itd_shiftice( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) )  ! Shift jl+1=>jl
629            ! Reset shift parameters
630            jdonor(1:npti,jl) = 0
631            zdaice(1:npti,jl) = 0._wp
632            zdvice(1:npti,jl) = 0._wp
633         ENDIF
634         !
635      END DO
636      !
637   END SUBROUTINE ice_itd_reb
638
639
640   SUBROUTINE ice_itd_init
641      !!------------------------------------------------------------------
642      !!                ***  ROUTINE ice_itd_init ***
643      !!
644      !! ** Purpose :   Initializes the ice thickness distribution
645      !! ** Method  :   ...
646      !! ** input   :   Namelist namitd
647      !!-------------------------------------------------------------------
648      INTEGER  ::   jl            ! dummy loop index
649      INTEGER  ::   ios, ioptio   ! Local integer output status for namelist read
650      REAL(wp) ::   zhmax, znum, zden, zalpha   !   -      -
651      !
652      NAMELIST/namitd/ ln_cat_hfn, rn_himean, ln_cat_usr, rn_catbnd, rn_himin
653      !!------------------------------------------------------------------
654      !
655      REWIND( numnam_ice_ref )      ! Namelist namitd in reference namelist : Parameters for ice
656      READ  ( numnam_ice_ref, namitd, IOSTAT = ios, ERR = 901)
657901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namitd in reference namelist', lwp )
658      REWIND( numnam_ice_cfg )      ! Namelist namitd in configuration namelist : Parameters for ice
659      READ  ( numnam_ice_cfg, namitd, IOSTAT = ios, ERR = 902 )
660902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namitd in configuration namelist', lwp )
661      IF(lwm) WRITE( numoni, namitd )
662      !
663      IF(lwp) THEN                  ! control print
664         WRITE(numout,*)
665         WRITE(numout,*) 'ice_itd_init: Initialization of ice cat distribution '
666         WRITE(numout,*) '~~~~~~~~~~~~'
667         WRITE(numout,*) '   Namelist namitd: '
668         WRITE(numout,*) '      Ice categories are defined by a function of rn_himean**(-0.05)    ln_cat_hfn = ', ln_cat_hfn
669         WRITE(numout,*) '         mean ice thickness in the domain                               rn_himean  = ', rn_himean
670         WRITE(numout,*) '      Ice categories are defined by rn_catbnd                           ln_cat_usr = ', ln_cat_usr
671         WRITE(numout,*) '      minimum ice thickness                                             rn_himin   = ', rn_himin 
672      ENDIF
673      !
674      !-----------------------------------!
675      !  Thickness categories boundaries  !
676      !-----------------------------------!
677      !                             !== set the choice of ice categories ==!
678      ioptio = 0 
679      IF( ln_cat_hfn ) THEN   ;   ioptio = ioptio + 1   ;   nice_catbnd = np_cathfn    ;   ENDIF
680      IF( ln_cat_usr ) THEN   ;   ioptio = ioptio + 1   ;   nice_catbnd = np_catusr    ;   ENDIF
681      IF( ioptio /= 1 )   CALL ctl_stop( 'ice_itd_init: choose one and only one ice categories boundaries' )
682      !
683      SELECT CASE( nice_catbnd )
684      !                                !------------------------!
685      CASE( np_cathfn )                ! h^(-alpha) function
686         !                             !------------------------!
687         zalpha = 0.05_wp
688         zhmax  = 3._wp * rn_himean
689         hi_max(0) = 0._wp
690         DO jl = 1, jpl
691            znum = jpl * ( zhmax+1 )**zalpha
692            zden = REAL( jpl-jl , wp ) * ( zhmax + 1._wp )**zalpha + REAL( jl , wp )
693            hi_max(jl) = ( znum / zden )**(1./zalpha) - 1
694         END DO
695         !                             !------------------------!
696      CASE( np_catusr )                ! user defined
697         !                             !------------------------!
698         DO jl = 0, jpl
699            hi_max(jl) = rn_catbnd(jl)
700         END DO
701         !
702      END SELECT
703      !
704      DO jl = 1, jpl                ! mean thickness by category
705         hi_mean(jl) = ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp
706      END DO
707      !
708      hi_max(jpl) = 99._wp          ! set to a big value to ensure that all ice is thinner than hi_max(jpl)
709      !
710      IF(lwp) WRITE(numout,*)
711      IF(lwp) WRITE(numout,*) '   ===>>>   resulting thickness category boundaries :'
712      IF(lwp) WRITE(numout,*) '            hi_max(:)= ', hi_max(0:jpl)
713      !
714      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')
715      !
716   END SUBROUTINE ice_itd_init
717
718#else
719   !!----------------------------------------------------------------------
720   !!   Default option :         Empty module         NO SI3 sea-ice model
721   !!----------------------------------------------------------------------
722#endif
723
724   !!======================================================================
725END MODULE iceitd
Note: See TracBrowser for help on using the repository browser.