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

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

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

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