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

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

changes in style - part6 - one more round

File size: 33.9 KB
Line 
1MODULE iceitd
2   !!======================================================================
3   !!                       ***  MODULE iceitd ***
4   !!   LIM3 ice model : ice thickness distribution: Thermodynamics
5   !!======================================================================
6   !! History :   -   !          (W. H. Lipscomb and E.C. Hunke) CICE (c) original code
7   !!            3.0  ! 2005-12  (M. Vancoppenolle) adaptation to LIM-3
8   !!             -   ! 2006-06  (M. Vancoppenolle) adaptation to include salt, age
9   !!             -   ! 2007-04  (M. Vancoppenolle) Mass conservation checked
10   !!----------------------------------------------------------------------
11#if defined key_lim3
12   !!----------------------------------------------------------------------
13   !!   'key_lim3'                                       LIM3 sea-ice model
14   !!----------------------------------------------------------------------
15   !!   ice_itd_rem   :
16   !!   ice_itd_reb   :
17   !!   ice_itd_glinear  :
18   !!   ice_itd_shiftice :
19   !!----------------------------------------------------------------------
20   USE par_oce        ! ocean parameters
21   USE dom_oce        ! ocean domain
22   USE phycst         ! physical constants
23   USE ice1D          ! sea-ice: thermodynamic variables
24   USE ice            ! sea-ice: variables
25   USE icectl         ! sea-ice: conservation tests
26   USE icetab         ! sea-ice: convert 1D<=>2D
27   !
28   USE prtctl         ! Print control
29   USE in_out_manager ! I/O manager
30   USE lib_mpp        ! MPP library
31   USE lib_fortran    ! to use key_nosignedzero
32
33   IMPLICIT NONE
34   PRIVATE
35
36   PUBLIC   ice_itd_init  ! called in icestp
37   PUBLIC   ice_itd_rem   ! called in icethd
38   PUBLIC   ice_itd_reb   ! called in iceerr
39
40   ! ** ice-thickness distribution namelist (namice_itd) **
41   REAL(wp) ::   rn_himean        ! mean thickness of the domain (used to compute the distribution)
42
43   !!----------------------------------------------------------------------
44   !! NEMO/ICE 4.0 , NEMO Consortium (2017)
45   !! $Id: iceitd.F90 8420 2017-08-08 12:18:46Z clem $
46   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
47   !!----------------------------------------------------------------------
48CONTAINS
49
50   SUBROUTINE ice_itd_rem( kt )
51      !!------------------------------------------------------------------
52      !!                ***  ROUTINE ice_itd_rem ***
53      !!
54      !! ** Purpose :   computes the redistribution of ice thickness
55      !!              after thermodynamic growth of ice thickness
56      !!
57      !! ** Method  : Linear remapping
58      !!
59      !! References : W.H. Lipscomb, JGR 2001
60      !!------------------------------------------------------------------
61      INTEGER , INTENT (in) ::   kt      ! Ocean time step
62      !
63      INTEGER  ::   ji, jj, jl, jcat     ! dummy loop index
64      INTEGER  ::   nidx2                ! local integer
65      REAL(wp) ::   zx1, zwk1, zdh0, zetamin, zdamax   ! local scalars
66      REAL(wp) ::   zx2, zwk2, zda0, zetamax           !   -      -
67      REAL(wp) ::   zx3       
68      REAL(wp) ::   zslope          ! used to compute local thermodynamic "speeds"
69
70      INTEGER , DIMENSION(jpij)       ::   idxice2         ! compute remapping or not
71      INTEGER , DIMENSION(jpij,jpl-1) ::   jdonor          ! donor category index
72      REAL(wp), DIMENSION(jpij,jpl)   ::   zdhice          ! ice thickness increment
73      REAL(wp), DIMENSION(jpij,jpl)   ::   g0, g1          ! coefficients for fitting the line of the ITD
74      REAL(wp), DIMENSION(jpij,jpl)   ::   hL, hR          ! left and right boundary for the ITD for each thickness
75      REAL(wp), DIMENSION(jpij,jpl-1) ::   zdaice, zdvice  ! local increment of ice area and volume
76      REAL(wp), DIMENSION(jpij)       ::   zhb0, zhb1      ! category boundaries for thinnes categories
77      REAL(wp), DIMENSION(jpij,0:jpl) ::   zhbnew          ! new boundaries of ice categories
78      !!------------------------------------------------------------------
79
80      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_itd_rem: remapping ice thickness distribution' 
81
82      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'iceitd_rem', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
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 ice_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 ice_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      !  4) Compute g(h)
184      !-----------------------------------------------------------------------------------------------
185      IF( nidx > 0 ) THEN
186         !
187         zhb0(:) = hi_max(0)   ;   zhb1(:) = hi_max(1)
188         g0(:,:) = 0._wp       ;   g1(:,:) = 0._wp 
189         hL(:,:) = 0._wp       ;   hR(:,:) = 0._wp 
190         !
191         DO jl = 1, jpl
192            !
193            CALL tab_2d_1d( nidx, idxice(1:nidx), ht_ib_1d(1:nidx), ht_i_b(:,:,jl) )
194            CALL tab_2d_1d( nidx, idxice(1:nidx), ht_i_1d (1:nidx), ht_i(:,:,jl)   )
195            CALL tab_2d_1d( nidx, idxice(1:nidx), a_i_1d  (1:nidx), a_i(:,:,jl)    )
196            CALL tab_2d_1d( nidx, idxice(1:nidx), v_i_1d  (1:nidx), v_i(:,:,jl)    )
197            !
198            IF( jl == 1 ) THEN
199               
200               ! --- g(h) for category 1 --- !
201               CALL ice_itd_glinear( zhb0(1:nidx)  , zhb1(1:nidx)  , ht_ib_1d(1:nidx)  , a_i_1d(1:nidx)  ,  &   ! in
202                  &                  g0  (1:nidx,1), g1  (1:nidx,1), hL      (1:nidx,1), hR    (1:nidx,1)   )   ! out
203                  !
204               ! Area lost due to melting of thin ice
205               DO ji = 1, nidx
206                  !
207                  IF( a_i_1d(ji) > epsi10 ) THEN
208                     !
209                     zdh0 =  ht_i_1d(ji) - ht_ib_1d(ji)               
210                     IF( zdh0 < 0.0 ) THEN      !remove area from category 1
211                        zdh0 = MIN( -zdh0, hi_max(1) )
212                        !Integrate g(1) from 0 to dh0 to estimate area melted
213                        zetamax = MIN( zdh0, hR(ji,1) ) - hL(ji,1)
214                        !
215                        IF( zetamax > 0.0 ) THEN
216                           zx1    = zetamax
217                           zx2    = 0.5 * zetamax * zetamax 
218                           zda0   = g1(ji,1) * zx2 + g0(ji,1) * zx1                        ! ice area removed
219                           zdamax = a_i_1d(ji) * (1.0 - ht_i_1d(ji) / ht_ib_1d(ji) ) ! Constrain new thickness <= ht_i               
220                           zda0   = MIN( zda0, zdamax )                                                  ! ice area lost due to melting
221                           !     of thin ice (zdamax > 0)
222                           ! Remove area, conserving volume
223                           ht_i_1d(ji) = ht_i_1d(ji) * a_i_1d(ji) / ( a_i_1d(ji) - zda0 )
224                           a_i_1d(ji)  = a_i_1d(ji) - zda0
225                           v_i_1d(ji)  = a_i_1d(ji) * ht_i_1d(ji) ! clem-useless ?
226                        ENDIF
227                        !
228                     ELSE ! if ice accretion zdh0 > 0
229                        ! zhbnew was 0, and is shifted to the right to account for thin ice growth in openwater (F0 = f1)
230                        zhbnew(ji,0) = MIN( zdh0, hi_max(1) ) 
231                     ENDIF
232                     !
233                  ENDIF
234                  !
235               END DO
236               !
237               CALL tab_1d_2d( nidx, idxice(1:nidx), ht_i_1d (1:nidx), ht_i(:,:,jl)   )
238               CALL tab_1d_2d( nidx, idxice(1:nidx), a_i_1d  (1:nidx), a_i(:,:,jl)    )
239               CALL tab_1d_2d( nidx, idxice(1:nidx), v_i_1d  (1:nidx), v_i(:,:,jl)    )
240               !
241            ENDIF ! jl=1
242            !
243            ! --- g(h) for each thickness category --- ! 
244            CALL ice_itd_glinear( zhbnew(1:nidx,jl-1), zhbnew(1:nidx,jl), ht_i_1d(1:nidx)   , a_i_1d(1:nidx)   ,  &   ! in
245               &                  g0    (1:nidx,jl  ), g1    (1:nidx,jl), hL     (1:nidx,jl), hR    (1:nidx,jl)   )   ! out
246            !
247         END DO
248         
249         !-----------------------------------------------------------------------------------------------
250         !  5) Compute area and volume to be shifted across each boundary (Eq. 18)
251         !-----------------------------------------------------------------------------------------------
252         DO jl = 1, jpl - 1
253            !
254            DO ji = 1, nidx
255               !
256               ! left and right integration limits in eta space
257               IF (zhbnew(ji,jl) > hi_max(jl)) THEN ! Hn* > Hn => transfer from jl to jl+1
258                  zetamin = MAX( hi_max(jl)   , hL(ji,jl) ) - hL(ji,jl)   ! hi_max(jl) - hL
259                  zetamax = MIN( zhbnew(ji,jl), hR(ji,jl) ) - hL(ji,jl)   ! hR - hL
260                  jdonor(ji,jl) = jl
261               ELSE                                 ! Hn* <= Hn => transfer from jl+1 to jl
262                  zetamin = 0.0
263                  zetamax = MIN( hi_max(jl), hR(ji,jl+1) ) - hL(ji,jl+1)  ! hi_max(jl) - hL
264                  jdonor(ji,jl) = jl + 1
265               ENDIF
266               zetamax = MAX( zetamax, zetamin ) ! no transfer if etamax < etamin
267               !
268               zx1  = zetamax - zetamin
269               zwk1 = zetamin * zetamin
270               zwk2 = zetamax * zetamax
271               zx2  = 0.5 * ( zwk2 - zwk1 )
272               zwk1 = zwk1 * zetamin
273               zwk2 = zwk2 * zetamax
274               zx3  = 1.0 / 3.0 * ( zwk2 - zwk1 )
275               jcat = jdonor(ji,jl)
276               zdaice(ji,jl) = g1(ji,jcat)*zx2 + g0(ji,jcat)*zx1
277               zdvice(ji,jl) = g1(ji,jcat)*zx3 + g0(ji,jcat)*zx2 + zdaice(ji,jl)*hL(ji,jcat)
278               !
279            END DO
280         END DO
281         
282         !----------------------------------------------------------------------------------------------
283         ! 6) Shift ice between categories
284         !----------------------------------------------------------------------------------------------
285         CALL ice_itd_shiftice ( jdonor(1:nidx,:), zdaice(1:nidx,:), zdvice(1:nidx,:) )
286         
287         !----------------------------------------------------------------------------------------------
288         ! 7) Make sure ht_i >= minimum ice thickness hi_min
289         !----------------------------------------------------------------------------------------------
290         CALL tab_2d_1d( nidx, idxice(1:nidx), ht_i_1d (1:nidx), ht_i(:,:,1)   )
291         CALL tab_2d_1d( nidx, idxice(1:nidx), a_i_1d  (1:nidx), a_i(:,:,1)    )
292         CALL tab_2d_1d( nidx, idxice(1:nidx), a_ip_1d (1:nidx), a_ip(:,:,1)   )
293         
294         DO ji = 1, nidx
295            IF ( a_i_1d(ji) > epsi10 .AND. ht_i_1d(ji) < rn_himin ) THEN
296               a_i_1d (ji) = a_i_1d(ji) * ht_i_1d(ji) / rn_himin 
297               ! MV MP 2016
298               IF ( nn_pnd_scheme > 0 ) THEN
299                  a_ip_1d(ji) = a_ip_1d(ji) * ht_i_1d(ji) / rn_himin
300               ENDIF
301               ! END MV MP 2016
302               ht_i_1d(ji) = rn_himin
303            ENDIF
304         END DO
305         !
306         CALL tab_1d_2d( nidx, idxice(1:nidx), ht_i_1d (1:nidx), ht_i(:,:,1)   )
307         CALL tab_1d_2d( nidx, idxice(1:nidx), a_i_1d  (1:nidx), a_i(:,:,1)    )
308         CALL tab_1d_2d( nidx, idxice(1:nidx), a_ip_1d (1:nidx), a_ip(:,:,1)   )
309         !
310      ENDIF
311      !
312      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'iceitd_rem', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
313      !
314   END SUBROUTINE ice_itd_rem
315
316
317   SUBROUTINE ice_itd_glinear( HbL, Hbr, phice, paice, pg0, pg1, phL, phR )
318      !!------------------------------------------------------------------
319      !!                ***  ROUTINE ice_itd_glinear ***
320      !!
321      !! ** Purpose :   build g(h) satisfying area and volume constraints (Eq. 6 and 9)
322      !!
323      !! ** Method  :   g(h) is linear and written as: g(eta) = g1(eta) + g0
324      !!                with eta = h - HL
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) ::   z1_3 , z2_3  ! 1/3 , 2/3
333      REAL(wp) ::   zh13         ! HbL + 1/3 * (HbR - HbL)
334      REAL(wp) ::   zh23         ! HbL + 2/3 * (HbR - HbL)
335      REAL(wp) ::   zdhr         ! 1 / (hR - hL)
336      REAL(wp) ::   zwk1, zwk2   ! temporary variables
337      !!------------------------------------------------------------------
338      !
339      z1_3 = 1._wp / 3._wp
340      z2_3 = 2._wp / 3._wp
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 = z1_3 * ( 2._wp * phL(ji) +         phR(ji) )
353            zh23 = z1_3 * (         phL(ji) + 2._wp * 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 * ( z2_3 - zwk2 )                    ! Eq. 14
364            pg1(ji) = 2._wp * zdhr * zwk1 * ( zwk2 - 0.5_wp )   ! 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 ice_itd_glinear
376
377
378   SUBROUTINE ice_itd_shiftice( kdonor, pdaice, pdvice )
379      !!------------------------------------------------------------------
380      !!                ***  ROUTINE ice_itd_shiftice ***
381      !!
382      !! ** Purpose :   shift ice across category boundaries, conserving everything
383      !!              ( area, volume, energy, age*vol, and mass of salt )
384      !!------------------------------------------------------------------
385      INTEGER , DIMENSION(:,:), INTENT(in) ::   kdonor   ! donor category index
386      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pdaice   ! ice area transferred across boundary
387      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pdvice   ! ice volume transferred across boundary
388      !
389      INTEGER  ::   ji, jj, jl, jk     ! dummy loop indices
390      INTEGER  ::   ii, ij, jl2, jl1   ! local integers
391      REAL(wp) ::   ztrans             ! ice/snow transferred
392      REAL(wp), DIMENSION(jpij)     ::   zworka, zworkv   ! workspace
393      REAL(wp), DIMENSION(jpij,jpl) ::   zaTsfn           !  -    -
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            !
423            IF( jl1 > 0 ) THEN
424               !
425               IF ( jl1 == jl  ) THEN   ;   jl2 = jl1+1
426               ELSE                     ;   jl2 = jl 
427               ENDIF
428               !
429               IF( v_i_2d(ji,jl1) >= epsi10 ) THEN   ;   zworkv(ji) = pdvice(ji,jl) / v_i_2d(ji,jl1)
430               ELSE                                  ;   zworkv(ji) = 0._wp
431               ENDIF
432               IF( a_i_2d(ji,jl1) >= epsi10 ) THEN   ;   zworka(ji) = pdaice(ji,jl) / a_i_2d(ji,jl1)
433               ELSE                                  ;   zworka(ji) = 0._wp
434               ENDIF
435               !
436               a_i_2d(ji,jl1) = a_i_2d(ji,jl1) - pdaice(ji,jl)       ! Ice areas
437               a_i_2d(ji,jl2) = a_i_2d(ji,jl2) + pdaice(ji,jl)
438               !
439               v_i_2d(ji,jl1) = v_i_2d(ji,jl1) - pdvice(ji,jl)       ! Ice volumes
440               v_i_2d(ji,jl2) = v_i_2d(ji,jl2) + pdvice(ji,jl)
441               !
442               ztrans         = v_s_2d(ji,jl1) * zworkv(ji)          ! Snow volumes
443               v_s_2d(ji,jl1) = v_s_2d(ji,jl1) - ztrans
444               v_s_2d(ji,jl2) = v_s_2d(ji,jl2) + ztrans 
445               !         
446               !                                                     ! Ice age
447               ztrans             = oa_i_2d(ji,jl1) * pdaice(ji,jl)  !!clem: should be * zworka(ji) but it does not work ????
448               oa_i_2d(ji,jl1)    = oa_i_2d(ji,jl1) - ztrans
449               oa_i_2d(ji,jl2)    = oa_i_2d(ji,jl2) + ztrans
450               !
451               ztrans             = smv_i_2d(ji,jl1) * zworkv(ji)    ! Ice salinity
452               !
453               smv_i_2d(ji,jl1)   = smv_i_2d(ji,jl1) - ztrans
454               smv_i_2d(ji,jl2)   = smv_i_2d(ji,jl2) + ztrans
455               !
456               !                                                     ! Surface temperature
457               ztrans          = t_su_2d(ji,jl1) * pdaice(ji,jl)     !!clem: should be * zworka(ji) but it does not work ????
458               zaTsfn(ji,jl1)  = zaTsfn(ji,jl1) - ztrans
459               zaTsfn(ji,jl2)  = zaTsfn(ji,jl2) + ztrans
460               
461               ! MV MP 2016
462               IF ( nn_pnd_scheme > 0 ) THEN
463                  !                                                  ! Pond fraction
464                  ztrans          = a_ip_2d(ji,jl1) * pdaice(ji,jl) !!clem: should be * zworka(ji) but it does not work
465                  a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - ztrans
466                  a_ip_2d(ji,jl2) = a_ip_2d(ji,jl2) + ztrans
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            ENDIF   ! jl1 >0
475         END DO
476         !
477         DO jk = 1, nlay_s         !--- Snow heat content
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               !
485               IF( jl1 > 0 ) THEN
486                  IF(jl1 == jl) THEN  ;  jl2 = jl+1
487                  ELSE                ;  jl2 = jl
488                  ENDIF
489                  !
490                  ztrans            = e_s(ii,ij,jk,jl1) * zworkv(ji)
491                  e_s(ii,ij,jk,jl1) = e_s(ii,ij,jk,jl1) - ztrans
492                  e_s(ii,ij,jk,jl2) = e_s(ii,ij,jk,jl2) + ztrans
493               ENDIF
494            END DO
495         END DO
496
497         DO jk = 1, nlay_i         !--- Ice heat content
498            DO ji = 1, nidx
499               ii = MOD( idxice(ji) - 1, jpi ) + 1
500               ij = ( idxice(ji) - 1 ) / jpi + 1
501               !
502               jl1 = kdonor(ji,jl)
503               !
504               IF( jl1 > 0 ) THEN
505                  IF(jl1 == jl) THEN  ;  jl2 = jl+1
506                  ELSE                ;  jl2 = jl
507                  ENDIF
508                  !
509                  ztrans            = e_i(ii,ij,jk,jl1) * zworkv(ji)
510                  e_i(ii,ij,jk,jl1) = e_i(ii,ij,jk,jl1) - ztrans
511                  e_i(ii,ij,jk,jl2) = e_i(ii,ij,jk,jl2) + ztrans
512               ENDIF
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      WHERE( a_i_2d(1:nidx,:) >= epsi20 )
522         ht_i_2d(1:nidx,:)  =  v_i_2d(1:nidx,:) / a_i_2d(1:nidx,:) 
523         t_su_2d(1:nidx,:)  =  zaTsfn(1:nidx,:) / a_i_2d(1:nidx,:) 
524      ELSEWHERE
525         ht_i_2d(1:nidx,:)  = 0._wp
526         t_su_2d(1:nidx,:)  = rt0
527      END WHERE
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   END SUBROUTINE ice_itd_shiftice
540   
541
542   SUBROUTINE ice_itd_reb( kt )
543      !!------------------------------------------------------------------
544      !!                ***  ROUTINE ice_itd_reb ***
545      !!
546      !! ** Purpose : rebin - rebins thicknesses into defined categories
547      !!
548      !! ** Method  : If a category thickness is out of bounds, shift part (for down to top)
549      !!              or entire (for top to down) area, volume, and energy
550      !!              to the neighboring category
551      !!------------------------------------------------------------------
552      INTEGER , INTENT (in) ::   kt      ! Ocean time step
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      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_itd_reb: rebining ice thickness distribution' 
560
561      jdonor(:,:) = 0
562      zdaice(:,:) = 0._wp
563      zdvice(:,:) = 0._wp
564      !
565      !                       !---------------------------------------
566      DO jl = 1, jpl-1        ! identify thicknesses that are too big
567         !                    !---------------------------------------
568         nidx = 0   ;   idxice(:) = 0
569         DO jj = 1, jpj
570            DO ji = 1, jpi
571               IF( a_i(ji,jj,jl) > epsi10 .AND. v_i(ji,jj,jl) > (a_i(ji,jj,jl) * hi_max(jl)) ) THEN
572                  nidx = nidx + 1
573                  idxice( nidx ) = (jj - 1) * jpi + ji                 
574               ENDIF
575            END DO
576         END DO
577         !
578!!clem   CALL tab_2d_1d( nidx, idxice(1:nidx), ht_i_1d (1:nidx), ht_i(:,:,jl)   )
579         CALL tab_2d_1d( nidx, idxice(1:nidx), a_i_1d  (1:nidx), a_i(:,:,jl)    )
580         CALL tab_2d_1d( nidx, idxice(1:nidx), v_i_1d  (1:nidx), v_i(:,:,jl)    )
581         !
582         DO ji = 1, nidx
583            jdonor(ji,jl)  = jl 
584            ! how much of a_i you send in cat sup is somewhat arbitrary
585!!clem: these do not work properly after a restart (I do not know why)
586!!          zdaice(ji,jl)  = a_i_1d(ji) * ( ht_i_1d(ji) - hi_max(jl) + epsi10 ) / ht_i_1d(ji) 
587!!          zdvice(ji,jl)  = v_i_1d(ji) - ( a_i_1d(ji) - zdaice(ji,jl) ) * ( hi_max(jl) - epsi10 )
588!!clem: these do not work properly after a restart (I do not know why)
589!!            zdaice(ji,jl)  = a_i_1d(ji)
590!!            zdvice(ji,jl)  = v_i_1d(ji)
591!!clem: these are from UCL and work ok
592            zdaice(ji,jl)  = a_i_1d(ji) * 0.5_wp
593            zdvice(ji,jl)  = v_i_1d(ji) - zdaice(ji,jl) * ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp
594         END DO
595         !
596         IF( nidx > 0 ) THEN
597            CALL ice_itd_shiftice( jdonor(1:nidx,:), zdaice(1:nidx,:), zdvice(1:nidx,:) )  ! Shift jl=>jl+1
598            ! Reset shift parameters
599            jdonor(1:nidx,jl) = 0
600            zdaice(1:nidx,jl) = 0._wp
601            zdvice(1:nidx,jl) = 0._wp
602         ENDIF
603         !
604      END DO
605
606      !                       !-----------------------------------------
607      DO jl = jpl-1, 1, -1    ! Identify thicknesses that are too small
608         !                    !-----------------------------------------
609         nidx = 0 ; idxice(:) = 0
610         DO jj = 1, jpj
611            DO ji = 1, jpi
612               IF( a_i(ji,jj,jl+1) > epsi10 .AND. v_i(ji,jj,jl+1) <= (a_i(ji,jj,jl+1) * hi_max(jl)) ) THEN
613                  nidx = nidx + 1
614                  idxice( nidx ) = (jj - 1) * jpi + ji                 
615               ENDIF
616            END DO
617         END DO
618         !
619         CALL tab_2d_1d( nidx, idxice(1:nidx), a_i_1d  (1:nidx), a_i(:,:,jl+1)    ) ! jl+1 is ok
620         CALL tab_2d_1d( nidx, idxice(1:nidx), v_i_1d  (1:nidx), v_i(:,:,jl+1)    ) ! jl+1 is ok
621         DO ji = 1, nidx
622            jdonor(ji,jl) = jl + 1
623            zdaice(ji,jl) = a_i_1d(ji) 
624            zdvice(ji,jl) = v_i_1d(ji)
625         END DO
626         !
627         IF( nidx > 0 ) THEN
628            CALL ice_itd_shiftice( jdonor(1:nidx,:), zdaice(1:nidx,:), zdvice(1:nidx,:) )  ! Shift jl+1=>jl
629            ! Reset shift parameters
630            jdonor(1:nidx,jl) = 0
631            zdaice(1:nidx,jl) = 0._wp
632            zdvice(1:nidx,jl) = 0._wp
633         ENDIF
634         !
635      END DO
636      !
637   END SUBROUTINE ice_itd_reb
638
639   SUBROUTINE ice_itd_init
640      !!------------------------------------------------------------------
641      !!                ***  ROUTINE ice_itd_init ***
642      !!
643      !! ** Purpose :   Initializes the ice thickness distribution
644      !! ** Method  :   ...
645      !! ** input   :   Namelist namice_itd
646      !!-------------------------------------------------------------------
647      INTEGER  ::   jl    ! dummy loop index
648      INTEGER  ::   ios   ! Local integer output status for namelist read
649      REAL(wp) ::   zhmax, znum, zden, zalpha   !   -      -
650      !!
651      NAMELIST/namice_itd/ rn_himean, rn_himin
652      !!------------------------------------------------------------------
653      !
654      REWIND( numnam_ice_ref )      ! Namelist namice_itd in reference namelist : Parameters for ice
655      READ  ( numnam_ice_ref, namice_itd, IOSTAT = ios, ERR = 901)
656901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namice_itd in reference namelist', lwp )
657
658      REWIND( numnam_ice_cfg )      ! Namelist namice_itd in configuration namelist : Parameters for ice
659      READ  ( numnam_ice_cfg, namice_itd, IOSTAT = ios, ERR = 902 )
660902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namice_itd in configuration namelist', lwp )
661      IF(lwm) WRITE ( numoni, namice_itd )
662      !
663      IF(lwp) THEN                  ! control print
664         WRITE(numout,*)
665         WRITE(numout,*) 'ice_itd_init: Initialization of ice cat distribution '
666         WRITE(numout,*) '~~~~~~~~~~~~'
667         WRITE(numout,*) '   Namelist namice_itd : '
668         WRITE(numout,*) '      mean ice thickness in the domain               rn_himean = ', rn_himean
669         WRITE(numout,*) '      minimum ice thickness                          rn_himin  = ', rn_himin 
670      ENDIF
671      !
672      !-----------------------------------!
673      !  Thickness categories boundaries  !
674      !-----------------------------------!
675      !
676      zalpha = 0.05_wp              ! max of each category (from h^(-alpha) function)
677      zhmax  = 3._wp * rn_himean
678      DO jl = 1, jpl
679         znum = jpl * ( zhmax+1 )**zalpha
680         zden = REAL( jpl-jl , wp ) * ( zhmax + 1._wp )**zalpha + REAL( jl , wp )
681         hi_max(jl) = ( znum / zden )**(1./zalpha) - 1
682      END DO
683      !
684      DO jl = 1, jpl                ! mean thickness by category
685         hi_mean(jl) = ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp
686      END DO
687      !
688      hi_max(jpl) = 99._wp          ! set to a big value to ensure that all ice is thinner than hi_max(jpl)
689      !
690      IF(lwp) WRITE(numout,*)
691      IF(lwp) WRITE(numout,*) '   ===>>>   resulting thickness category boundaries :'
692      IF(lwp) WRITE(numout,*) '            hi_max(:)= ', hi_max(0:jpl)
693      !
694   END SUBROUTINE ice_itd_init
695
696#else
697   !!----------------------------------------------------------------------
698   !!   Default option :         Empty module          NO LIM sea-ice model
699   !!----------------------------------------------------------------------
700#endif
701
702   !!======================================================================
703END MODULE iceitd
Note: See TracBrowser for help on using the repository browser.