source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90 @ 8410

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

correction on the last commit

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