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 @ 8426

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

last routine names to be changed
File size: 29.9 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 (ocean directory)
23   USE ice1D            ! LIM-3 thermodynamic variables
24   USE ice              ! LIM-3 variables
25   USE icectl           ! conservation tests
26   USE icetab
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/LIM3 4.0 , UCL - NEMO Consortium (2010)
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      !!------------------------------------------------------------------
329      REAL(wp), DIMENSION(:), INTENT(in   ) ::   HbL, HbR      ! left and right category boundaries
330      REAL(wp), DIMENSION(:), INTENT(in   ) ::   phice, paice  ! ice thickness and concentration
331      REAL(wp), DIMENSION(:), INTENT(inout) ::   pg0, pg1      ! coefficients in linear equation for g(eta)
332      REAL(wp), DIMENSION(:), INTENT(inout) ::   phL, phR      ! min and max value of range over which g(h) > 0
333      !
334      INTEGER  ::   ji           ! horizontal indices
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         DO ji = 1, nidx
342            !
343            IF( paice(ji) > epsi10  .AND. phice(ji) > 0._wp )  THEN
344
345               ! Initialize hL and hR
346               phL(ji) = HbL(ji)
347               phR(ji) = HbR(ji)
348
349               ! Change hL or hR if hice falls outside central third of range,
350               ! so that hice is in the central third of the range [HL HR]
351               zh13 = 1.0 / 3.0 * ( 2.0 * phL(ji) +       phR(ji) )
352               zh23 = 1.0 / 3.0 * (       phL(ji) + 2.0 * phR(ji) )
353
354               IF    ( phice(ji) < zh13 ) THEN   ;   phR(ji) = 3._wp * phice(ji) - 2._wp * phL(ji) ! move HR to the left
355               ELSEIF( phice(ji) > zh23 ) THEN   ;   phL(ji) = 3._wp * phice(ji) - 2._wp * phR(ji) ! move HL to the right
356               ENDIF
357
358               ! Compute coefficients of g(eta) = g0 + g1*eta
359               zdhr = 1._wp / (phR(ji) - phL(ji))
360               zwk1 = 6._wp * paice(ji) * zdhr
361               zwk2 = ( phice(ji) - phL(ji) ) * zdhr
362               pg0(ji) = zwk1 * ( 2._wp / 3._wp - zwk2 )      ! Eq. 14
363               pg1(ji) = 2._wp * zdhr * zwk1 * ( zwk2 - 0.5 ) ! Eq. 14
364               !
365            ELSE  ! remap_flag = .false. or a_i < epsi10
366               phL(ji) = 0._wp
367               phR(ji) = 0._wp
368               pg0(ji) = 0._wp
369               pg1(ji) = 0._wp
370            ENDIF
371            !
372         END DO
373      !
374   END SUBROUTINE ice_itd_glinear
375
376
377   SUBROUTINE ice_itd_shiftice( kdonor, pdaice, pdvice )
378      !!------------------------------------------------------------------
379      !!                ***  ROUTINE ice_itd_shiftice ***
380      !!
381      !! ** Purpose :   shift ice across category boundaries, conserving everything
382      !!              ( area, volume, energy, age*vol, and mass of salt )
383      !!
384      !! ** Method  :
385      !!------------------------------------------------------------------
386      INTEGER , DIMENSION(:,:), INTENT(in) ::   kdonor   ! donor category index
387      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pdaice   ! ice area transferred across boundary
388      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pdvice   ! ice volume transferred across boundary
389
390      INTEGER ::   ii,ij, ji, jj, jl, jl2, jl1, jk   ! dummy loop indices
391
392      REAL(wp), DIMENSION(jpij,jpl) ::   zaTsfn
393      REAL(wp), DIMENSION(jpij)     ::   zworka            ! temporary array used here
394      REAL(wp), DIMENSION(jpij)     ::   zworkv            ! temporary array used here
395
396      REAL(wp) ::   ztrans             ! ice/snow transferred
397      !!------------------------------------------------------------------
398         
399      CALL tab_3d_2d( nidx, idxice(1:nidx), ht_i_2d (1:nidx,1:jpl), ht_i   )
400      CALL tab_3d_2d( nidx, idxice(1:nidx), a_i_2d  (1:nidx,1:jpl), a_i    )
401      CALL tab_3d_2d( nidx, idxice(1:nidx), v_i_2d  (1:nidx,1:jpl), v_i    )
402      CALL tab_3d_2d( nidx, idxice(1:nidx), v_s_2d  (1:nidx,1:jpl), v_s    )
403      CALL tab_3d_2d( nidx, idxice(1:nidx), oa_i_2d (1:nidx,1:jpl), oa_i   )
404      CALL tab_3d_2d( nidx, idxice(1:nidx), smv_i_2d(1:nidx,1:jpl), smv_i  )
405      CALL tab_3d_2d( nidx, idxice(1:nidx), a_ip_2d (1:nidx,1:jpl), a_ip   )
406      CALL tab_3d_2d( nidx, idxice(1:nidx), v_ip_2d (1:nidx,1:jpl), v_ip   )
407      CALL tab_3d_2d( nidx, idxice(1:nidx), t_su_2d (1:nidx,1:jpl), t_su   )
408
409      !----------------------------------------------------------------------------------------------
410      ! 1) Define a variable equal to a_i*T_su
411      !----------------------------------------------------------------------------------------------
412      DO jl = 1, jpl
413         DO ji = 1, nidx
414            zaTsfn(ji,jl) = a_i_2d(ji,jl) * t_su_2d(ji,jl)
415         END DO
416      END DO
417     
418      !-------------------------------------------------------------------------------
419      ! 2) Transfer volume and energy between categories
420      !-------------------------------------------------------------------------------
421      DO jl = 1, jpl - 1
422         DO ji = 1, nidx
423           
424            jl1 = kdonor(ji,jl)
425            IF    ( jl1 == jl  ) THEN   ;   jl2 = jl1+1
426            ELSE                        ;   jl2 = jl 
427            ENDIF
428
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            ! Ice areas
436            a_i_2d(ji,jl1) = a_i_2d(ji,jl1) - pdaice(ji,jl)
437            a_i_2d(ji,jl2) = a_i_2d(ji,jl2) + pdaice(ji,jl)
438               
439            ! Ice volumes
440            v_i_2d(ji,jl1) = v_i_2d(ji,jl1) - pdvice(ji,jl) 
441            v_i_2d(ji,jl2) = v_i_2d(ji,jl2) + pdvice(ji,jl)
442           
443            ! Snow volumes
444            ztrans         = v_s_2d(ji,jl1) * zworkv(ji)
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            ! Ice salinity
454            ztrans             = smv_i_2d(ji,jl1) * zworkv(ji)
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               
470               ! Pond volume (also proportional to da/a)
471               ztrans          = v_ip_2d(ji,jl1) * pdaice(ji,jl) !!clem: should be * zworka(ji) but it does not work
472               v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - ztrans
473               v_ip_2d(ji,jl2) = v_ip_2d(ji,jl2) + ztrans
474            ENDIF
475            ! END MV MP 2016
476               
477         END DO
478
479         ! Snow heat content
480         DO jk = 1, nlay_s
481
482            DO ji = 1, nidx
483               ii = MOD( idxice(ji) - 1, jpi ) + 1
484               ij = ( idxice(ji) - 1 ) / jpi + 1
485
486               jl1 = kdonor(ji,jl)
487               IF(jl1 == jl) THEN  ;  jl2 = jl+1
488               ELSE                ;  jl2 = jl
489               ENDIF
490               
491               ztrans            = e_s(ii,ij,jk,jl1) * zworkv(ji)
492               e_s(ii,ij,jk,jl1) = e_s(ii,ij,jk,jl1) - ztrans
493               e_s(ii,ij,jk,jl2) = e_s(ii,ij,jk,jl2) + ztrans 
494            END DO
495         END DO
496
497         
498         ! Ice heat content
499         DO jk = 1, nlay_i
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      !
543   END SUBROUTINE ice_itd_shiftice
544   
545
546   SUBROUTINE ice_itd_reb
547      !!------------------------------------------------------------------
548      !!                ***  ROUTINE ice_itd_reb ***
549      !!
550      !! ** Purpose : rebin - rebins thicknesses into defined categories
551      !!
552      !! ** Method  : If a category thickness is out of bounds, shift part (for down to top)
553      !!              or entire (for top to down) area, volume, and energy
554      !!              to the neighboring category
555      !!------------------------------------------------------------------
556      INTEGER ::   ji, jj, jl   ! dummy loop indices
557      !
558      INTEGER , DIMENSION(jpij,jpl-1) ::   jdonor           ! donor category index
559      REAL(wp), DIMENSION(jpij,jpl-1) ::   zdaice, zdvice   ! ice area and volume transferred
560      !!------------------------------------------------------------------
561     
562      DO jl = 1, jpl-1
563         jdonor(:,jl) = 0
564         zdaice(:,jl) = 0._wp
565         zdvice(:,jl) = 0._wp
566      END DO
567
568      !---------------------------------------
569      ! identify thicknesses that are too big
570      !---------------------------------------
571      DO jl = 1, jpl - 1
572
573         nidx = 0 ; idxice(:) = 0
574         DO jj = 1, jpj
575            DO ji = 1, jpi
576               IF( a_i(ji,jj,jl) > epsi10 .AND. v_i(ji,jj,jl) > (a_i(ji,jj,jl) * hi_max(jl)) ) THEN
577                  nidx = nidx + 1
578                  idxice( nidx ) = (jj - 1) * jpi + ji                 
579               ENDIF
580            ENDDO
581         ENDDO
582         
583!!clem   CALL tab_2d_1d( nidx, idxice(1:nidx), ht_i_1d (1:nidx), ht_i(:,:,jl)   )
584         CALL tab_2d_1d( nidx, idxice(1:nidx), a_i_1d  (1:nidx), a_i(:,:,jl)    )
585         CALL tab_2d_1d( nidx, idxice(1:nidx), v_i_1d  (1:nidx), v_i(:,:,jl)    )
586
587         DO ji = 1, nidx
588            jdonor(ji,jl)  = jl 
589            ! how much of a_i you send in cat sup is somewhat arbitrary
590!!clem: these do not work properly after a restart (I do not know why)
591!!          zdaice(ji,jl)  = a_i_1d(ji) * ( ht_i_1d(ji) - hi_max(jl) + epsi10 ) / ht_i_1d(ji) 
592!!          zdvice(ji,jl)  = v_i_1d(ji) - ( a_i_1d(ji) - zdaice(ji,jl) ) * ( hi_max(jl) - epsi10 )
593!!clem: these do not work properly after a restart (I do not know why)
594!!            zdaice(ji,jl)  = a_i_1d(ji)
595!!            zdvice(ji,jl)  = v_i_1d(ji)
596!!clem: these are from UCL and work ok
597            zdaice(ji,jl)  = a_i_1d(ji) * 0.5_wp
598            zdvice(ji,jl)  = v_i_1d(ji) - zdaice(ji,jl) * ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp
599           
600         END DO
601
602         IF( nidx > 0 ) THEN
603            CALL ice_itd_shiftice( jdonor(1:nidx,:), zdaice(1:nidx,:), zdvice(1:nidx,:) )  ! Shift jl=>jl+1
604            ! Reset shift parameters
605            jdonor(1:nidx,jl) = 0
606            zdaice(1:nidx,jl) = 0._wp
607            zdvice(1:nidx,jl) = 0._wp
608         ENDIF
609         !
610      END DO
611
612      !-----------------------------------------
613      ! Identify thicknesses that are too small
614      !-----------------------------------------
615      DO jl = jpl - 1, 1, -1
616
617         nidx = 0 ; idxice(:) = 0
618         DO jj = 1, jpj
619            DO ji = 1, jpi
620               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
621                  nidx = nidx + 1
622                  idxice( nidx ) = (jj - 1) * jpi + ji                 
623               ENDIF
624            ENDDO
625         ENDDO
626
627         CALL tab_2d_1d( nidx, idxice(1:nidx), a_i_1d  (1:nidx), a_i(:,:,jl+1)    ) ! jl+1 is ok
628         CALL tab_2d_1d( nidx, idxice(1:nidx), v_i_1d  (1:nidx), v_i(:,:,jl+1)    ) ! jl+1 is ok
629         DO ji = 1, nidx
630            jdonor(ji,jl) = jl + 1
631            zdaice(ji,jl) = a_i_1d(ji) 
632            zdvice(ji,jl) = v_i_1d(ji)
633         END DO
634         
635         IF( nidx > 0 ) THEN
636            CALL ice_itd_shiftice( jdonor(1:nidx,:), zdaice(1:nidx,:), zdvice(1:nidx,:) )  ! Shift jl+1=>jl
637            ! Reset shift parameters
638            jdonor(1:nidx,jl) = 0
639            zdaice(1:nidx,jl) = 0._wp
640            zdvice(1:nidx,jl) = 0._wp
641         ENDIF
642
643      END DO
644      !
645   END SUBROUTINE ice_itd_reb
646
647#endif
648   !!======================================================================
649END MODULE iceitd
Note: See TracBrowser for help on using the repository browser.