source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90 @ 8514

Last change on this file since 8514 was 8514, checked in by clem, 3 years ago

changes in style - part5 - almost done

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