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 branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

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

Last change on this file since 8498 was 8498, checked in by clem, 7 years ago

changes in style - part2 -
File size: 30.7 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            IF ( jl1 == jl  ) THEN   ;   jl2 = jl1+1
425            ELSE                     ;   jl2 = jl 
426            ENDIF
427            !
428            IF( v_i_2d(ji,jl1) >= epsi10 ) THEN   ;   zworkv(ji) = pdvice(ji,jl) / v_i_2d(ji,jl1)
429            ELSE                                  ;   zworkv(ji) = 0._wp
430            ENDIF
431            IF( a_i_2d(ji,jl1) >= epsi10 ) THEN   ;   zworka(ji) = pdaice(ji,jl) / a_i_2d(ji,jl1)
432            ELSE                                  ;   zworka(ji) = 0._wp
433            ENDIF
434            !
435            a_i_2d(ji,jl1) = a_i_2d(ji,jl1) - pdaice(ji,jl)       ! Ice areas
436            a_i_2d(ji,jl2) = a_i_2d(ji,jl2) + pdaice(ji,jl)
437            !
438            v_i_2d(ji,jl1) = v_i_2d(ji,jl1) - pdvice(ji,jl)       ! Ice volumes
439            v_i_2d(ji,jl2) = v_i_2d(ji,jl2) + pdvice(ji,jl)
440            !
441            ztrans         = v_s_2d(ji,jl1) * zworkv(ji)          ! Snow volumes
442            v_s_2d(ji,jl1) = v_s_2d(ji,jl1) - ztrans
443            v_s_2d(ji,jl2) = v_s_2d(ji,jl2) + ztrans 
444            !         
445            !                                                     ! Ice age
446            ztrans             = oa_i_2d(ji,jl1) * pdaice(ji,jl)  !!clem: should be * zworka(ji) but it does not work ????
447            oa_i_2d(ji,jl1)    = oa_i_2d(ji,jl1) - ztrans
448            oa_i_2d(ji,jl2)    = oa_i_2d(ji,jl2) + ztrans
449            !
450            ztrans             = smv_i_2d(ji,jl1) * zworkv(ji)    ! Ice salinity
451
452            smv_i_2d(ji,jl1)   = smv_i_2d(ji,jl1) - ztrans
453            smv_i_2d(ji,jl2)   = smv_i_2d(ji,jl2) + ztrans
454            !
455            !                                                     ! Surface temperature
456            ztrans          = t_su_2d(ji,jl1) * pdaice(ji,jl)     !!clem: should be * zworka(ji) but it does not work ????
457            zaTsfn(ji,jl1)  = zaTsfn(ji,jl1) - ztrans
458            zaTsfn(ji,jl2)  = zaTsfn(ji,jl2) + ztrans
459            ! 
460            ! MV MP 2016
461            IF ( nn_pnd_scheme > 0 ) THEN
462               !                                                  ! Pond fraction
463               ztrans          = a_ip_2d(ji,jl1) * pdaice(ji,jl) !!clem: should be * zworka(ji) but it does not work
464               a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - ztrans
465               a_ip_2d(ji,jl2) = a_ip_2d(ji,jl2) + ztrans
466               !                                                  ! Pond volume (also proportional to da/a)
467               ztrans          = v_ip_2d(ji,jl1) * pdaice(ji,jl) !!clem: should be * zworka(ji) but it does not work
468               v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - ztrans
469               v_ip_2d(ji,jl2) = v_ip_2d(ji,jl2) + ztrans
470            ENDIF
471            ! END MV MP 2016
472         END DO
473         !
474         DO jk = 1, nlay_s         !--- Snow heat content
475            !
476            DO ji = 1, nidx
477               ii = MOD( idxice(ji) - 1, jpi ) + 1
478               ij = ( idxice(ji) - 1 ) / jpi + 1
479               !
480               jl1 = kdonor(ji,jl)
481               IF(jl1 == jl) THEN  ;  jl2 = jl+1
482               ELSE                ;  jl2 = jl
483               ENDIF
484               !
485               ztrans            = e_s(ii,ij,jk,jl1) * zworkv(ji)
486               e_s(ii,ij,jk,jl1) = e_s(ii,ij,jk,jl1) - ztrans
487               e_s(ii,ij,jk,jl2) = e_s(ii,ij,jk,jl2) + ztrans 
488            END DO
489         END DO
490
491         DO jk = 1, nlay_i         !--- Ice heat content
492            DO ji = 1, nidx
493               ii = MOD( idxice(ji) - 1, jpi ) + 1
494               ij = ( idxice(ji) - 1 ) / jpi + 1
495               !
496               jl1 = kdonor(ji,jl)
497               IF(jl1 == jl) THEN  ;  jl2 = jl+1
498               ELSE                ;  jl2 = jl
499               ENDIF
500               !
501               ztrans            = e_i(ii,ij,jk,jl1) * zworkv(ji)
502               e_i(ii,ij,jk,jl1) = e_i(ii,ij,jk,jl1) - ztrans
503               e_i(ii,ij,jk,jl2) = e_i(ii,ij,jk,jl2) + ztrans 
504            END DO
505         END DO
506         !
507      END DO                   ! boundaries, 1 to jpl-1
508     
509      !-------------------------------------------------------------------------------
510      ! 3) Update ice thickness and temperature
511      !-------------------------------------------------------------------------------
512      DO jl = 1, jpl
513         DO ji = 1, nidx
514            IF ( a_i_2d(ji,jl) > epsi10 ) THEN
515               ht_i_2d(ji,jl)  =  v_i_2d(ji,jl) / a_i_2d(ji,jl) 
516               t_su_2d(ji,jl)  =  zaTsfn(ji,jl) / a_i_2d(ji,jl) 
517            ELSE
518               ht_i_2d(ji,jl)  = 0._wp
519               t_su_2d(ji,jl)  = rt0
520            ENDIF
521         END DO
522      END DO
523      !
524      CALL tab_2d_3d( nidx, idxice(1:nidx), ht_i_2d (1:nidx,1:jpl), ht_i  )
525      CALL tab_2d_3d( nidx, idxice(1:nidx), a_i_2d  (1:nidx,1:jpl), a_i   )
526      CALL tab_2d_3d( nidx, idxice(1:nidx), v_i_2d  (1:nidx,1:jpl), v_i   )
527      CALL tab_2d_3d( nidx, idxice(1:nidx), v_s_2d  (1:nidx,1:jpl), v_s   )
528      CALL tab_2d_3d( nidx, idxice(1:nidx), oa_i_2d (1:nidx,1:jpl), oa_i  )
529      CALL tab_2d_3d( nidx, idxice(1:nidx), smv_i_2d(1:nidx,1:jpl), smv_i )
530      CALL tab_2d_3d( nidx, idxice(1:nidx), a_ip_2d (1:nidx,1:jpl), a_ip  )
531      CALL tab_2d_3d( nidx, idxice(1:nidx), v_ip_2d (1:nidx,1:jpl), v_ip  )
532      CALL tab_2d_3d( nidx, idxice(1:nidx), t_su_2d (1:nidx,1:jpl), t_su  )
533      !
534   END SUBROUTINE ice_itd_shiftice
535   
536
537   SUBROUTINE ice_itd_reb
538      !!------------------------------------------------------------------
539      !!                ***  ROUTINE ice_itd_reb ***
540      !!
541      !! ** Purpose : rebin - rebins thicknesses into defined categories
542      !!
543      !! ** Method  : If a category thickness is out of bounds, shift part (for down to top)
544      !!              or entire (for top to down) area, volume, and energy
545      !!              to the neighboring category
546      !!------------------------------------------------------------------
547      INTEGER ::   ji, jj, jl   ! dummy loop indices
548      !
549      INTEGER , DIMENSION(jpij,jpl-1) ::   jdonor           ! donor category index
550      REAL(wp), DIMENSION(jpij,jpl-1) ::   zdaice, zdvice   ! ice area and volume transferred
551      !!------------------------------------------------------------------
552      !
553      jdonor(:,:) = 0
554      zdaice(:,:) = 0._wp
555      zdvice(:,:) = 0._wp
556      !
557      !                       !---------------------------------------
558      DO jl = 1, jpl-1        ! identify thicknesses that are too big
559         !                    !---------------------------------------
560         nidx = 0   ;   idxice(:) = 0
561         DO jj = 1, jpj
562            DO ji = 1, jpi
563               IF( a_i(ji,jj,jl) > epsi10 .AND. v_i(ji,jj,jl) > (a_i(ji,jj,jl) * hi_max(jl)) ) THEN
564                  nidx = nidx + 1
565                  idxice( nidx ) = (jj - 1) * jpi + ji                 
566               ENDIF
567            END DO
568         END DO
569         !
570!!clem   CALL tab_2d_1d( nidx, idxice(1:nidx), ht_i_1d (1:nidx), ht_i(:,:,jl)   )
571         CALL tab_2d_1d( nidx, idxice(1:nidx), a_i_1d  (1:nidx), a_i(:,:,jl)    )
572         CALL tab_2d_1d( nidx, idxice(1:nidx), v_i_1d  (1:nidx), v_i(:,:,jl)    )
573         !
574         DO ji = 1, nidx
575            jdonor(ji,jl)  = jl 
576            ! how much of a_i you send in cat sup is somewhat arbitrary
577!!clem: these do not work properly after a restart (I do not know why)
578!!          zdaice(ji,jl)  = a_i_1d(ji) * ( ht_i_1d(ji) - hi_max(jl) + epsi10 ) / ht_i_1d(ji) 
579!!          zdvice(ji,jl)  = v_i_1d(ji) - ( a_i_1d(ji) - zdaice(ji,jl) ) * ( hi_max(jl) - epsi10 )
580!!clem: these do not work properly after a restart (I do not know why)
581!!            zdaice(ji,jl)  = a_i_1d(ji)
582!!            zdvice(ji,jl)  = v_i_1d(ji)
583!!clem: these are from UCL and work ok
584            zdaice(ji,jl)  = a_i_1d(ji) * 0.5_wp
585            zdvice(ji,jl)  = v_i_1d(ji) - zdaice(ji,jl) * ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp
586         END DO
587         !
588         IF( nidx > 0 ) THEN
589            CALL ice_itd_shiftice( jdonor(1:nidx,:), zdaice(1:nidx,:), zdvice(1:nidx,:) )  ! Shift jl=>jl+1
590            ! Reset shift parameters
591            jdonor(1:nidx,jl) = 0
592            zdaice(1:nidx,jl) = 0._wp
593            zdvice(1:nidx,jl) = 0._wp
594         ENDIF
595         !
596      END DO
597
598      !                       !-----------------------------------------
599      DO jl = jpl-1, 1, -1    ! Identify thicknesses that are too small
600         !                    !-----------------------------------------
601         nidx = 0 ; idxice(:) = 0
602         DO jj = 1, jpj
603            DO ji = 1, jpi
604               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
605                  nidx = nidx + 1
606                  idxice( nidx ) = (jj - 1) * jpi + ji                 
607               ENDIF
608            END DO
609         END DO
610         !
611         CALL tab_2d_1d( nidx, idxice(1:nidx), a_i_1d  (1:nidx), a_i(:,:,jl+1)    ) ! jl+1 is ok
612         CALL tab_2d_1d( nidx, idxice(1:nidx), v_i_1d  (1:nidx), v_i(:,:,jl+1)    ) ! jl+1 is ok
613         DO ji = 1, nidx
614            jdonor(ji,jl) = jl + 1
615            zdaice(ji,jl) = a_i_1d(ji) 
616            zdvice(ji,jl) = v_i_1d(ji)
617         END DO
618         !
619         IF( nidx > 0 ) THEN
620            CALL ice_itd_shiftice( jdonor(1:nidx,:), zdaice(1:nidx,:), zdvice(1:nidx,:) )  ! Shift jl+1=>jl
621            ! Reset shift parameters
622            jdonor(1:nidx,jl) = 0
623            zdaice(1:nidx,jl) = 0._wp
624            zdvice(1:nidx,jl) = 0._wp
625         ENDIF
626         !
627      END DO
628      !
629   END SUBROUTINE ice_itd_reb
630
631#else
632   !!----------------------------------------------------------------------
633   !!   Default option :         Empty module          NO LIM sea-ice model
634   !!----------------------------------------------------------------------
635#endif
636
637   !!======================================================================
638END MODULE iceitd
Note: See TracBrowser for help on using the repository browser.