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

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

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