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

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

changes in style - part1 - (now the code looks better txs to Gurvan's comments)

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