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

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

changes in style - part3 - move output into proper files and correct a bug in debug mode

File size: 31.0 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_rem   ! called in icethd
37   PUBLIC   ice_itd_reb   ! called in iceerr
38
39   !!----------------------------------------------------------------------
40   !! NEMO/ICE 4.0 , NEMO Consortium (2017)
41   !! $Id: iceitd.F90 8420 2017-08-08 12:18:46Z clem $
42   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
43   !!----------------------------------------------------------------------
44CONTAINS
45
46   SUBROUTINE ice_itd_rem( kt )
47      !!------------------------------------------------------------------
48      !!                ***  ROUTINE ice_itd_rem ***
49      !!
50      !! ** Purpose :   computes the redistribution of ice thickness
51      !!              after thermodynamic growth of ice thickness
52      !!
53      !! ** Method  : Linear remapping
54      !!
55      !! References : W.H. Lipscomb, JGR 2001
56      !!------------------------------------------------------------------
57      INTEGER , INTENT (in) ::   kt      ! Ocean time step
58      !
59      INTEGER  ::   ji, jj, jl, jcat     ! dummy loop index
60      INTEGER  ::   nidx2                ! local integer
61      REAL(wp) ::   zx1, zwk1, zdh0, zetamin, zdamax   ! local scalars
62      REAL(wp) ::   zx2, zwk2, zda0, zetamax           !   -      -
63      REAL(wp) ::   zx3       
64      REAL(wp) ::   zslope          ! used to compute local thermodynamic "speeds"
65      REAL(wp) ::   zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b ! conservation check
66
67      INTEGER , DIMENSION(jpij)       ::   idxice2         ! compute remapping or not
68      INTEGER , DIMENSION(jpij,jpl-1) ::   jdonor          ! donor category index
69      REAL(wp), DIMENSION(jpij,jpl)   ::   zdhice          ! ice thickness increment
70      REAL(wp), DIMENSION(jpij,jpl)   ::   g0, g1          ! coefficients for fitting the line of the ITD
71      REAL(wp), DIMENSION(jpij,jpl)   ::   hL, hR          ! left and right boundary for the ITD for each thickness
72      REAL(wp), DIMENSION(jpij,jpl-1) ::   zdaice, zdvice  ! local increment of ice area and volume
73      REAL(wp), DIMENSION(jpij)       ::   zhb0, zhb1      ! category boundaries for thinnes categories
74      REAL(wp), DIMENSION(jpij,0:jpl) ::   zhbnew          ! new boundaries of ice categories
75      !!------------------------------------------------------------------
76
77      IF( kt == nit000 .AND. lwp) THEN
78         WRITE(numout,*)
79         WRITE(numout,*) 'ice_itd_rem  : Remapping the ice thickness distribution'
80         WRITE(numout,*) '~~~~~~~~~~~~~~~'
81      ENDIF
82
83      IF( ln_limdiachk ) 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      !-----------------------------------------------------------------------------------------------
185      !  4) Compute g(h)
186      !-----------------------------------------------------------------------------------------------
187      IF( nidx > 0 ) THEN
188         !
189         zhb0(:) = hi_max(0)   ;   zhb1(:) = hi_max(1)
190         g0(:,:) = 0._wp       ;   g1(:,:) = 0._wp 
191         hL(:,:) = 0._wp       ;   hR(:,:) = 0._wp 
192         !
193         DO jl = 1, jpl
194            !
195            CALL tab_2d_1d( nidx, idxice(1:nidx), ht_ib_1d(1:nidx), ht_i_b(:,:,jl) )
196            CALL tab_2d_1d( nidx, idxice(1:nidx), ht_i_1d (1:nidx), ht_i(:,:,jl)   )
197            CALL tab_2d_1d( nidx, idxice(1:nidx), a_i_1d  (1:nidx), a_i(:,:,jl)    )
198            CALL tab_2d_1d( nidx, idxice(1:nidx), v_i_1d  (1:nidx), v_i(:,:,jl)    )
199            !
200            IF( jl == 1 ) THEN
201               
202               ! --- g(h) for category 1 --- !
203               CALL ice_itd_glinear( zhb0(1:nidx)  , zhb1(1:nidx)  , ht_ib_1d(1:nidx)  , a_i_1d(1:nidx)  ,  &   ! in
204                  &                  g0  (1:nidx,1), g1  (1:nidx,1), hL      (1:nidx,1), hR    (1:nidx,1)   )   ! out
205                  !
206               ! Area lost due to melting of thin ice
207               DO ji = 1, nidx
208                  !
209                  IF( a_i_1d(ji) > epsi10 ) THEN
210                     !
211                     zdh0 =  ht_i_1d(ji) - ht_ib_1d(ji)               
212                     IF( zdh0 < 0.0 ) THEN      !remove area from category 1
213                        zdh0 = MIN( -zdh0, hi_max(1) )
214                        !Integrate g(1) from 0 to dh0 to estimate area melted
215                        zetamax = MIN( zdh0, hR(ji,1) ) - hL(ji,1)
216                        !
217                        IF( zetamax > 0.0 ) THEN
218                           zx1    = zetamax
219                           zx2    = 0.5 * zetamax * zetamax 
220                           zda0   = g1(ji,1) * zx2 + g0(ji,1) * zx1                        ! ice area removed
221                           zdamax = a_i_1d(ji) * (1.0 - ht_i_1d(ji) / ht_ib_1d(ji) ) ! Constrain new thickness <= ht_i               
222                           zda0   = MIN( zda0, zdamax )                                                  ! ice area lost due to melting
223                           !     of thin ice (zdamax > 0)
224                           ! Remove area, conserving volume
225                           ht_i_1d(ji) = ht_i_1d(ji) * a_i_1d(ji) / ( a_i_1d(ji) - zda0 )
226                           a_i_1d(ji)  = a_i_1d(ji) - zda0
227                           v_i_1d(ji)  = a_i_1d(ji) * ht_i_1d(ji) ! clem-useless ?
228                        ENDIF
229                        !
230                     ELSE ! if ice accretion zdh0 > 0
231                        ! zhbnew was 0, and is shifted to the right to account for thin ice growth in openwater (F0 = f1)
232                        zhbnew(ji,0) = MIN( zdh0, hi_max(1) ) 
233                     ENDIF
234                     !
235                  ENDIF
236                  !
237               END DO
238               !
239               CALL tab_1d_2d( nidx, idxice(1:nidx), ht_i_1d (1:nidx), ht_i(:,:,jl)   )
240               CALL tab_1d_2d( nidx, idxice(1:nidx), a_i_1d  (1:nidx), a_i(:,:,jl)    )
241               CALL tab_1d_2d( nidx, idxice(1:nidx), v_i_1d  (1:nidx), v_i(:,:,jl)    )
242               !
243            ENDIF ! jl=1
244            !
245            ! --- g(h) for each thickness category --- ! 
246            CALL ice_itd_glinear( zhbnew(1:nidx,jl-1), zhbnew(1:nidx,jl), ht_i_1d(1:nidx)   , a_i_1d(1:nidx)   ,  &   ! in
247               &                  g0    (1:nidx,jl  ), g1    (1:nidx,jl), hL     (1:nidx,jl), hR    (1:nidx,jl)   )   ! out
248            !
249         END DO
250         
251         !-----------------------------------------------------------------------------------------------
252         !  5) Compute area and volume to be shifted across each boundary (Eq. 18)
253         !-----------------------------------------------------------------------------------------------
254         DO jl = 1, jpl - 1
255            !
256            DO ji = 1, nidx
257               !
258               ! left and right integration limits in eta space
259               IF (zhbnew(ji,jl) > hi_max(jl)) THEN ! Hn* > Hn => transfer from jl to jl+1
260                  zetamin = MAX( hi_max(jl)   , hL(ji,jl) ) - hL(ji,jl)   ! hi_max(jl) - hL
261                  zetamax = MIN( zhbnew(ji,jl), hR(ji,jl) ) - hL(ji,jl)   ! hR - hL
262                  jdonor(ji,jl) = jl
263               ELSE                                 ! Hn* <= Hn => transfer from jl+1 to jl
264                  zetamin = 0.0
265                  zetamax = MIN( hi_max(jl), hR(ji,jl+1) ) - hL(ji,jl+1)  ! hi_max(jl) - hL
266                  jdonor(ji,jl) = jl + 1
267               ENDIF
268               zetamax = MAX( zetamax, zetamin ) ! no transfer if etamax < etamin
269               !
270               zx1  = zetamax - zetamin
271               zwk1 = zetamin * zetamin
272               zwk2 = zetamax * zetamax
273               zx2  = 0.5 * ( zwk2 - zwk1 )
274               zwk1 = zwk1 * zetamin
275               zwk2 = zwk2 * zetamax
276               zx3  = 1.0 / 3.0 * ( zwk2 - zwk1 )
277               jcat = jdonor(ji,jl)
278               zdaice(ji,jl) = g1(ji,jcat)*zx2 + g0(ji,jcat)*zx1
279               zdvice(ji,jl) = g1(ji,jcat)*zx3 + g0(ji,jcat)*zx2 + zdaice(ji,jl)*hL(ji,jcat)
280               !
281            END DO
282         END DO
283         
284         !----------------------------------------------------------------------------------------------
285         ! 6) Shift ice between categories
286         !----------------------------------------------------------------------------------------------
287         CALL ice_itd_shiftice ( jdonor(1:nidx,:), zdaice(1:nidx,:), zdvice(1:nidx,:) )
288         
289         !----------------------------------------------------------------------------------------------
290         ! 7) Make sure ht_i >= minimum ice thickness hi_min
291         !----------------------------------------------------------------------------------------------
292         CALL tab_2d_1d( nidx, idxice(1:nidx), ht_i_1d (1:nidx), ht_i(:,:,1)   )
293         CALL tab_2d_1d( nidx, idxice(1:nidx), a_i_1d  (1:nidx), a_i(:,:,1)    )
294         CALL tab_2d_1d( nidx, idxice(1:nidx), a_ip_1d (1:nidx), a_ip(:,:,1)    )
295         
296         DO ji = 1, nidx
297            IF ( a_i_1d(ji) > epsi10 .AND. ht_i_1d(ji) < rn_himin ) THEN
298               a_i_1d (ji) = a_i_1d(ji) * ht_i_1d(ji) / rn_himin 
299               ! MV MP 2016
300               IF ( nn_pnd_scheme > 0 ) THEN
301                  a_ip_1d(ji) = a_ip_1d(ji) * ht_i_1d(ji) / rn_himin
302               ENDIF
303               ! END MV MP 2016
304               ht_i_1d(ji) = rn_himin
305            ENDIF
306         END DO
307         !
308         CALL tab_1d_2d( nidx, idxice(1:nidx), ht_i_1d (1:nidx), ht_i(:,:,1)   )
309         CALL tab_1d_2d( nidx, idxice(1:nidx), a_i_1d  (1:nidx), a_i(:,:,1)    )
310         CALL tab_1d_2d( nidx, idxice(1:nidx), a_ip_1d (1:nidx), a_ip(:,:,1)    )
311         !
312      ENDIF
313      !
314      IF( ln_limdiachk )   CALL ice_cons_hsm(1, 'iceitd_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
315      !
316   END SUBROUTINE ice_itd_rem
317
318
319   SUBROUTINE ice_itd_glinear( HbL, Hbr, phice, paice, pg0, pg1, phL, phR )
320      !!------------------------------------------------------------------
321      !!                ***  ROUTINE ice_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, nidx
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 ice_itd_glinear
378
379
380   SUBROUTINE ice_itd_shiftice( kdonor, pdaice, pdvice )
381      !!------------------------------------------------------------------
382      !!                ***  ROUTINE ice_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( nidx, idxice(1:nidx), ht_i_2d (1:nidx,1:jpl), ht_i   )
399      CALL tab_3d_2d( nidx, idxice(1:nidx), a_i_2d  (1:nidx,1:jpl), a_i    )
400      CALL tab_3d_2d( nidx, idxice(1:nidx), v_i_2d  (1:nidx,1:jpl), v_i    )
401      CALL tab_3d_2d( nidx, idxice(1:nidx), v_s_2d  (1:nidx,1:jpl), v_s    )
402      CALL tab_3d_2d( nidx, idxice(1:nidx), oa_i_2d (1:nidx,1:jpl), oa_i   )
403      CALL tab_3d_2d( nidx, idxice(1:nidx), smv_i_2d(1:nidx,1:jpl), smv_i  )
404      CALL tab_3d_2d( nidx, idxice(1:nidx), a_ip_2d (1:nidx,1:jpl), a_ip   )
405      CALL tab_3d_2d( nidx, idxice(1:nidx), v_ip_2d (1:nidx,1:jpl), v_ip   )
406      CALL tab_3d_2d( nidx, idxice(1:nidx), t_su_2d (1:nidx,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, nidx
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, nidx
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) * pdaice(ji,jl)  !!clem: should be * zworka(ji) but it does not work ????
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             = smv_i_2d(ji,jl1) * zworkv(ji)    ! Ice salinity
454               !
455               smv_i_2d(ji,jl1)   = smv_i_2d(ji,jl1) - ztrans
456               smv_i_2d(ji,jl2)   = smv_i_2d(ji,jl2) + ztrans
457               !
458               !                                                     ! Surface temperature
459               ztrans          = t_su_2d(ji,jl1) * pdaice(ji,jl)     !!clem: should be * zworka(ji) but it does not work ????
460               zaTsfn(ji,jl1)  = zaTsfn(ji,jl1) - ztrans
461               zaTsfn(ji,jl2)  = zaTsfn(ji,jl2) + ztrans
462               
463               ! MV MP 2016
464               IF ( nn_pnd_scheme > 0 ) THEN
465                  !                                                  ! Pond fraction
466                  ztrans          = a_ip_2d(ji,jl1) * pdaice(ji,jl) !!clem: should be * zworka(ji) but it does not work
467                  a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - ztrans
468                  a_ip_2d(ji,jl2) = a_ip_2d(ji,jl2) + ztrans
469                  !                                                  ! Pond volume (also proportional to da/a)
470                  ztrans          = v_ip_2d(ji,jl1) * pdaice(ji,jl) !!clem: should be * zworka(ji) but it does not work
471                  v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - ztrans
472                  v_ip_2d(ji,jl2) = v_ip_2d(ji,jl2) + ztrans
473               ENDIF
474               ! END MV MP 2016
475               !
476            ENDIF   ! jl1 >0
477         END DO
478         !
479         DO jk = 1, nlay_s         !--- Snow heat content
480            !
481            DO ji = 1, nidx
482               ii = MOD( idxice(ji) - 1, jpi ) + 1
483               ij = ( idxice(ji) - 1 ) / jpi + 1
484               !
485               jl1 = kdonor(ji,jl)
486               !
487               IF( jl1 > 0 ) THEN
488                  IF(jl1 == jl) THEN  ;  jl2 = jl+1
489                  ELSE                ;  jl2 = jl
490                  ENDIF
491                  !
492                  ztrans            = e_s(ii,ij,jk,jl1) * zworkv(ji)
493                  e_s(ii,ij,jk,jl1) = e_s(ii,ij,jk,jl1) - ztrans
494                  e_s(ii,ij,jk,jl2) = e_s(ii,ij,jk,jl2) + ztrans
495               ENDIF
496            END DO
497         END DO
498
499         DO jk = 1, nlay_i         !--- Ice heat content
500            DO ji = 1, nidx
501               ii = MOD( idxice(ji) - 1, jpi ) + 1
502               ij = ( idxice(ji) - 1 ) / jpi + 1
503               !
504               jl1 = kdonor(ji,jl)
505               !
506               IF( jl1 > 0 ) THEN
507                  IF(jl1 == jl) THEN  ;  jl2 = jl+1
508                  ELSE                ;  jl2 = jl
509                  ENDIF
510                  !
511                  ztrans            = e_i(ii,ij,jk,jl1) * zworkv(ji)
512                  e_i(ii,ij,jk,jl1) = e_i(ii,ij,jk,jl1) - ztrans
513                  e_i(ii,ij,jk,jl2) = e_i(ii,ij,jk,jl2) + ztrans
514               ENDIF
515            END DO
516         END DO
517         !
518      END DO                   ! boundaries, 1 to jpl-1
519     
520      !-------------------------------------------------------------------------------
521      ! 3) Update ice thickness and temperature
522      !-------------------------------------------------------------------------------
523      WHERE( a_i_2d(1:nidx,:) >= epsi20 )
524         ht_i_2d(1:nidx,:)  =  v_i_2d(1:nidx,:) / a_i_2d(1:nidx,:) 
525         t_su_2d(1:nidx,:)  =  zaTsfn(1:nidx,:) / a_i_2d(1:nidx,:) 
526      ELSEWHERE
527         ht_i_2d(1:nidx,:)  = 0._wp
528         t_su_2d(1:nidx,:)  = rt0
529      END WHERE
530      !
531      CALL tab_2d_3d( nidx, idxice(1:nidx), ht_i_2d (1:nidx,1:jpl), ht_i  )
532      CALL tab_2d_3d( nidx, idxice(1:nidx), a_i_2d  (1:nidx,1:jpl), a_i   )
533      CALL tab_2d_3d( nidx, idxice(1:nidx), v_i_2d  (1:nidx,1:jpl), v_i   )
534      CALL tab_2d_3d( nidx, idxice(1:nidx), v_s_2d  (1:nidx,1:jpl), v_s   )
535      CALL tab_2d_3d( nidx, idxice(1:nidx), oa_i_2d (1:nidx,1:jpl), oa_i  )
536      CALL tab_2d_3d( nidx, idxice(1:nidx), smv_i_2d(1:nidx,1:jpl), smv_i )
537      CALL tab_2d_3d( nidx, idxice(1:nidx), a_ip_2d (1:nidx,1:jpl), a_ip  )
538      CALL tab_2d_3d( nidx, idxice(1:nidx), v_ip_2d (1:nidx,1:jpl), v_ip  )
539      CALL tab_2d_3d( nidx, idxice(1:nidx), t_su_2d (1:nidx,1:jpl), t_su  )
540      !
541   END SUBROUTINE ice_itd_shiftice
542   
543
544   SUBROUTINE ice_itd_reb
545      !!------------------------------------------------------------------
546      !!                ***  ROUTINE ice_itd_reb ***
547      !!
548      !! ** Purpose : rebin - rebins thicknesses into defined categories
549      !!
550      !! ** Method  : If a category thickness is out of bounds, shift part (for down to top)
551      !!              or entire (for top to down) area, volume, and energy
552      !!              to the neighboring category
553      !!------------------------------------------------------------------
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      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         nidx = 0   ;   idxice(:) = 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                  nidx = nidx + 1
572                  idxice( nidx ) = (jj - 1) * jpi + ji                 
573               ENDIF
574            END DO
575         END DO
576         !
577!!clem   CALL tab_2d_1d( nidx, idxice(1:nidx), ht_i_1d (1:nidx), ht_i(:,:,jl)   )
578         CALL tab_2d_1d( nidx, idxice(1:nidx), a_i_1d  (1:nidx), a_i(:,:,jl)    )
579         CALL tab_2d_1d( nidx, idxice(1:nidx), v_i_1d  (1:nidx), v_i(:,:,jl)    )
580         !
581         DO ji = 1, nidx
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) * ( ht_i_1d(ji) - hi_max(jl) + epsi10 ) / ht_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( nidx > 0 ) THEN
596            CALL ice_itd_shiftice( jdonor(1:nidx,:), zdaice(1:nidx,:), zdvice(1:nidx,:) )  ! Shift jl=>jl+1
597            ! Reset shift parameters
598            jdonor(1:nidx,jl) = 0
599            zdaice(1:nidx,jl) = 0._wp
600            zdvice(1:nidx,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         nidx = 0 ; idxice(:) = 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                  nidx = nidx + 1
613                  idxice( nidx ) = (jj - 1) * jpi + ji                 
614               ENDIF
615            END DO
616         END DO
617         !
618         CALL tab_2d_1d( nidx, idxice(1:nidx), a_i_1d  (1:nidx), a_i(:,:,jl+1)    ) ! jl+1 is ok
619         CALL tab_2d_1d( nidx, idxice(1:nidx), v_i_1d  (1:nidx), v_i(:,:,jl+1)    ) ! jl+1 is ok
620         DO ji = 1, nidx
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( nidx > 0 ) THEN
627            CALL ice_itd_shiftice( jdonor(1:nidx,:), zdaice(1:nidx,:), zdvice(1:nidx,:) )  ! Shift jl+1=>jl
628            ! Reset shift parameters
629            jdonor(1:nidx,jl) = 0
630            zdaice(1:nidx,jl) = 0._wp
631            zdvice(1:nidx,jl) = 0._wp
632         ENDIF
633         !
634      END DO
635      !
636   END SUBROUTINE ice_itd_reb
637
638#else
639   !!----------------------------------------------------------------------
640   !!   Default option :         Empty module          NO LIM sea-ice model
641   !!----------------------------------------------------------------------
642#endif
643
644   !!======================================================================
645END MODULE iceitd
Note: See TracBrowser for help on using the repository browser.