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

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

changes in style - part5 - start changing init routines

File size: 33.7 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 (namiceitd) **
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      REAL(wp) ::   zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b ! conservation check
70
71      INTEGER , DIMENSION(jpij)       ::   idxice2         ! compute remapping or not
72      INTEGER , DIMENSION(jpij,jpl-1) ::   jdonor          ! donor category index
73      REAL(wp), DIMENSION(jpij,jpl)   ::   zdhice          ! ice thickness increment
74      REAL(wp), DIMENSION(jpij,jpl)   ::   g0, g1          ! coefficients for fitting the line of the ITD
75      REAL(wp), DIMENSION(jpij,jpl)   ::   hL, hR          ! left and right boundary for the ITD for each thickness
76      REAL(wp), DIMENSION(jpij,jpl-1) ::   zdaice, zdvice  ! local increment of ice area and volume
77      REAL(wp), DIMENSION(jpij)       ::   zhb0, zhb1      ! category boundaries for thinnes categories
78      REAL(wp), DIMENSION(jpij,0:jpl) ::   zhbnew          ! new boundaries of ice categories
79      !!------------------------------------------------------------------
80
81      IF( kt == nit000 .AND. lwp) THEN
82         WRITE(numout,*)
83         WRITE(numout,*) 'ice_itd_rem  : Remapping the ice thickness distribution'
84         WRITE(numout,*) '~~~~~~~~~~~~~~~'
85      ENDIF
86
87      IF( ln_limdiachk ) CALL ice_cons_hsm(0, 'iceitd_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
88
89      !-----------------------------------------------------------------------------------------------
90      !  1) Identify grid cells with ice
91      !-----------------------------------------------------------------------------------------------
92      nidx = 0   ;   idxice(:) = 0
93      DO jj = 1, jpj
94         DO ji = 1, jpi
95            IF ( at_i(ji,jj) > epsi10 ) THEN
96               nidx         = nidx + 1
97               idxice( nidx ) = (jj - 1) * jpi + ji
98            ENDIF
99         END DO
100      END DO
101     
102      !-----------------------------------------------------------------------------------------------
103      !  2) Compute new category boundaries
104      !-----------------------------------------------------------------------------------------------
105      IF( nidx > 0 ) THEN
106
107         zdhice(:,:) = 0._wp
108         zhbnew(:,:) = 0._wp
109
110         CALL tab_3d_2d( nidx, idxice(1:nidx), ht_i_2d (1:nidx,1:jpl), ht_i   )
111         CALL tab_3d_2d( nidx, idxice(1:nidx), ht_ib_2d(1:nidx,1:jpl), ht_i_b )
112         CALL tab_3d_2d( nidx, idxice(1:nidx), a_i_2d  (1:nidx,1:jpl), a_i    )
113         CALL tab_3d_2d( nidx, idxice(1:nidx), a_ib_2d (1:nidx,1:jpl), a_i_b  )
114
115         DO jl = 1, jpl
116            ! Compute thickness change in each ice category
117            DO ji = 1, nidx
118               zdhice(ji,jl) = ht_i_2d(ji,jl) - ht_ib_2d(ji,jl)
119            END DO
120         END DO
121         
122         ! --- New boundaries for category 1:jpl-1 --- !
123         DO jl = 1, jpl - 1
124            !
125            DO ji = 1, nidx
126               !
127               ! --- New boundary: Hn* = Hn + Fn*dt --- !
128               !     Fn*dt = ( fn + (fn+1 - fn)/(hn+1 - hn) * (Hn - hn) ) * dt = zdhice + zslope * (Hmax - ht_i_b)
129               !
130               IF    ( a_ib_2d(ji,jl) >  epsi10 .AND. a_ib_2d(ji,jl+1) >  epsi10 ) THEN   ! a(jl+1) & a(jl) /= 0
131                  zslope           = ( zdhice(ji,jl+1) - zdhice(ji,jl) ) / ( ht_ib_2d(ji,jl+1) - ht_ib_2d(ji,jl) )
132                  zhbnew(ji,jl) = hi_max(jl) + zdhice(ji,jl) + zslope * ( hi_max(jl) - ht_ib_2d(ji,jl) )
133               ELSEIF( a_ib_2d(ji,jl) >  epsi10 .AND. a_ib_2d(ji,jl+1) <= epsi10 ) THEN   ! a(jl+1)=0 => Hn* = Hn + fn*dt
134                  zhbnew(ji,jl) = hi_max(jl) + zdhice(ji,jl)
135               ELSEIF( a_ib_2d(ji,jl) <= epsi10 .AND. a_ib_2d(ji,jl+1) >  epsi10 ) THEN   ! a(jl)=0 => Hn* = Hn + fn+1*dt
136                  zhbnew(ji,jl) = hi_max(jl) + zdhice(ji,jl+1)
137               ELSE                                                                       ! a(jl+1) & a(jl) = 0
138                  zhbnew(ji,jl) = hi_max(jl)
139               ENDIF
140               !
141               ! --- 2 conditions for remapping --- !
142               ! 1) hn(t+1)+espi < Hn* < hn+1(t+1)-epsi               
143               !    Note: hn(t+1) must not be too close to either HR or HL otherwise a division by nearly 0 is possible
144               !          in ice_itd_glinear in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice)
145               IF( a_i_2d(ji,jl  ) > epsi10 .AND. ht_i_2d(ji,jl  ) > ( zhbnew(ji,jl) - epsi10 ) )   idxice(ji) = 0
146               IF( a_i_2d(ji,jl+1) > epsi10 .AND. ht_i_2d(ji,jl+1) < ( zhbnew(ji,jl) + epsi10 ) )   idxice(ji) = 0
147               
148               ! 2) Hn-1 < Hn* < Hn+1 
149               IF( zhbnew(ji,jl) < hi_max(jl-1) )   idxice(ji) = 0
150               IF( zhbnew(ji,jl) > hi_max(jl+1) )   idxice(ji) = 0
151               
152            END DO
153         END DO
154         !
155         ! --- New boundaries for category jpl --- !
156         DO ji = 1, nidx
157            IF( a_i_2d(ji,jpl) > epsi10 ) THEN
158               zhbnew(ji,jpl) = MAX( hi_max(jpl-1), 3._wp * ht_i_2d(ji,jpl) - 2._wp * zhbnew(ji,jpl-1) )
159            ELSE
160               zhbnew(ji,jpl) = hi_max(jpl) 
161            ENDIF
162           
163            ! --- 1 additional condition for remapping (1st category) --- !
164            ! H0+epsi < h1(t) < H1-epsi
165            !    h1(t) must not be too close to either HR or HL otherwise a division by nearly 0 is possible
166            !    in ice_itd_glinear in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice)
167            IF( ht_ib_2d(ji,1) < ( hi_max(0) + epsi10 ) )   idxice(ji) = 0
168            IF( ht_ib_2d(ji,1) > ( hi_max(1) - epsi10 ) )   idxice(ji) = 0
169         END DO
170         !
171         !-----------------------------------------------------------------------------------------------
172         !  3) Identify cells where remapping
173         !-----------------------------------------------------------------------------------------------
174         nidx2 = 0 ; idxice2(:) = 0
175         DO ji = 1, nidx
176            IF( idxice(ji) /= 0 ) THEN
177               nidx2 = nidx2 + 1
178               idxice2(nidx2) = idxice(ji)
179               zhbnew(nidx2,:) = zhbnew(ji,:)  ! adjust zhbnew to new indices
180            ENDIF
181         END DO
182         idxice(:) = idxice2(:)
183         nidx      = nidx2
184         !
185      ENDIF
186   
187      !-----------------------------------------------------------------------------------------------
188      !  4) Compute g(h)
189      !-----------------------------------------------------------------------------------------------
190      IF( nidx > 0 ) THEN
191         !
192         zhb0(:) = hi_max(0)   ;   zhb1(:) = hi_max(1)
193         g0(:,:) = 0._wp       ;   g1(:,:) = 0._wp 
194         hL(:,:) = 0._wp       ;   hR(:,:) = 0._wp 
195         !
196         DO jl = 1, jpl
197            !
198            CALL tab_2d_1d( nidx, idxice(1:nidx), ht_ib_1d(1:nidx), ht_i_b(:,:,jl) )
199            CALL tab_2d_1d( nidx, idxice(1:nidx), ht_i_1d (1:nidx), ht_i(:,:,jl)   )
200            CALL tab_2d_1d( nidx, idxice(1:nidx), a_i_1d  (1:nidx), a_i(:,:,jl)    )
201            CALL tab_2d_1d( nidx, idxice(1:nidx), v_i_1d  (1:nidx), v_i(:,:,jl)    )
202            !
203            IF( jl == 1 ) THEN
204               
205               ! --- g(h) for category 1 --- !
206               CALL ice_itd_glinear( zhb0(1:nidx)  , zhb1(1:nidx)  , ht_ib_1d(1:nidx)  , a_i_1d(1:nidx)  ,  &   ! in
207                  &                  g0  (1:nidx,1), g1  (1:nidx,1), hL      (1:nidx,1), hR    (1:nidx,1)   )   ! out
208                  !
209               ! Area lost due to melting of thin ice
210               DO ji = 1, nidx
211                  !
212                  IF( a_i_1d(ji) > epsi10 ) THEN
213                     !
214                     zdh0 =  ht_i_1d(ji) - ht_ib_1d(ji)               
215                     IF( zdh0 < 0.0 ) THEN      !remove area from category 1
216                        zdh0 = MIN( -zdh0, hi_max(1) )
217                        !Integrate g(1) from 0 to dh0 to estimate area melted
218                        zetamax = MIN( zdh0, hR(ji,1) ) - hL(ji,1)
219                        !
220                        IF( zetamax > 0.0 ) THEN
221                           zx1    = zetamax
222                           zx2    = 0.5 * zetamax * zetamax 
223                           zda0   = g1(ji,1) * zx2 + g0(ji,1) * zx1                        ! ice area removed
224                           zdamax = a_i_1d(ji) * (1.0 - ht_i_1d(ji) / ht_ib_1d(ji) ) ! Constrain new thickness <= ht_i               
225                           zda0   = MIN( zda0, zdamax )                                                  ! ice area lost due to melting
226                           !     of thin ice (zdamax > 0)
227                           ! Remove area, conserving volume
228                           ht_i_1d(ji) = ht_i_1d(ji) * a_i_1d(ji) / ( a_i_1d(ji) - zda0 )
229                           a_i_1d(ji)  = a_i_1d(ji) - zda0
230                           v_i_1d(ji)  = a_i_1d(ji) * ht_i_1d(ji) ! clem-useless ?
231                        ENDIF
232                        !
233                     ELSE ! if ice accretion zdh0 > 0
234                        ! zhbnew was 0, and is shifted to the right to account for thin ice growth in openwater (F0 = f1)
235                        zhbnew(ji,0) = MIN( zdh0, hi_max(1) ) 
236                     ENDIF
237                     !
238                  ENDIF
239                  !
240               END DO
241               !
242               CALL tab_1d_2d( nidx, idxice(1:nidx), ht_i_1d (1:nidx), ht_i(:,:,jl)   )
243               CALL tab_1d_2d( nidx, idxice(1:nidx), a_i_1d  (1:nidx), a_i(:,:,jl)    )
244               CALL tab_1d_2d( nidx, idxice(1:nidx), v_i_1d  (1:nidx), v_i(:,:,jl)    )
245               !
246            ENDIF ! jl=1
247            !
248            ! --- g(h) for each thickness category --- ! 
249            CALL ice_itd_glinear( zhbnew(1:nidx,jl-1), zhbnew(1:nidx,jl), ht_i_1d(1:nidx)   , a_i_1d(1:nidx)   ,  &   ! in
250               &                  g0    (1:nidx,jl  ), g1    (1:nidx,jl), hL     (1:nidx,jl), hR    (1:nidx,jl)   )   ! out
251            !
252         END DO
253         
254         !-----------------------------------------------------------------------------------------------
255         !  5) Compute area and volume to be shifted across each boundary (Eq. 18)
256         !-----------------------------------------------------------------------------------------------
257         DO jl = 1, jpl - 1
258            !
259            DO ji = 1, nidx
260               !
261               ! left and right integration limits in eta space
262               IF (zhbnew(ji,jl) > hi_max(jl)) THEN ! Hn* > Hn => transfer from jl to jl+1
263                  zetamin = MAX( hi_max(jl)   , hL(ji,jl) ) - hL(ji,jl)   ! hi_max(jl) - hL
264                  zetamax = MIN( zhbnew(ji,jl), hR(ji,jl) ) - hL(ji,jl)   ! hR - hL
265                  jdonor(ji,jl) = jl
266               ELSE                                 ! Hn* <= Hn => transfer from jl+1 to jl
267                  zetamin = 0.0
268                  zetamax = MIN( hi_max(jl), hR(ji,jl+1) ) - hL(ji,jl+1)  ! hi_max(jl) - hL
269                  jdonor(ji,jl) = jl + 1
270               ENDIF
271               zetamax = MAX( zetamax, zetamin ) ! no transfer if etamax < etamin
272               !
273               zx1  = zetamax - zetamin
274               zwk1 = zetamin * zetamin
275               zwk2 = zetamax * zetamax
276               zx2  = 0.5 * ( zwk2 - zwk1 )
277               zwk1 = zwk1 * zetamin
278               zwk2 = zwk2 * zetamax
279               zx3  = 1.0 / 3.0 * ( zwk2 - zwk1 )
280               jcat = jdonor(ji,jl)
281               zdaice(ji,jl) = g1(ji,jcat)*zx2 + g0(ji,jcat)*zx1
282               zdvice(ji,jl) = g1(ji,jcat)*zx3 + g0(ji,jcat)*zx2 + zdaice(ji,jl)*hL(ji,jcat)
283               !
284            END DO
285         END DO
286         
287         !----------------------------------------------------------------------------------------------
288         ! 6) Shift ice between categories
289         !----------------------------------------------------------------------------------------------
290         CALL ice_itd_shiftice ( jdonor(1:nidx,:), zdaice(1:nidx,:), zdvice(1:nidx,:) )
291         
292         !----------------------------------------------------------------------------------------------
293         ! 7) Make sure ht_i >= minimum ice thickness hi_min
294         !----------------------------------------------------------------------------------------------
295         CALL tab_2d_1d( nidx, idxice(1:nidx), ht_i_1d (1:nidx), ht_i(:,:,1)   )
296         CALL tab_2d_1d( nidx, idxice(1:nidx), a_i_1d  (1:nidx), a_i(:,:,1)    )
297         CALL tab_2d_1d( nidx, idxice(1:nidx), a_ip_1d (1:nidx), a_ip(:,:,1)    )
298         
299         DO ji = 1, nidx
300            IF ( a_i_1d(ji) > epsi10 .AND. ht_i_1d(ji) < rn_himin ) THEN
301               a_i_1d (ji) = a_i_1d(ji) * ht_i_1d(ji) / rn_himin 
302               ! MV MP 2016
303               IF ( nn_pnd_scheme > 0 ) THEN
304                  a_ip_1d(ji) = a_ip_1d(ji) * ht_i_1d(ji) / rn_himin
305               ENDIF
306               ! END MV MP 2016
307               ht_i_1d(ji) = rn_himin
308            ENDIF
309         END DO
310         !
311         CALL tab_1d_2d( nidx, idxice(1:nidx), ht_i_1d (1:nidx), ht_i(:,:,1)   )
312         CALL tab_1d_2d( nidx, idxice(1:nidx), a_i_1d  (1:nidx), a_i(:,:,1)    )
313         CALL tab_1d_2d( nidx, idxice(1:nidx), a_ip_1d (1:nidx), a_ip(:,:,1)    )
314         !
315      ENDIF
316      !
317      IF( ln_limdiachk )   CALL ice_cons_hsm(1, 'iceitd_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
318      !
319   END SUBROUTINE ice_itd_rem
320
321
322   SUBROUTINE ice_itd_glinear( HbL, Hbr, phice, paice, pg0, pg1, phL, phR )
323      !!------------------------------------------------------------------
324      !!                ***  ROUTINE ice_itd_glinear ***
325      !!
326      !! ** Purpose :   build g(h) satisfying area and volume constraints (Eq. 6 and 9)
327      !!
328      !! ** Method  :   g(h) is linear and written as: g(eta) = g1(eta) + g0
329      !!                with eta = h - HL
330      !!------------------------------------------------------------------
331      REAL(wp), DIMENSION(:), INTENT(in   ) ::   HbL, HbR      ! left and right category boundaries
332      REAL(wp), DIMENSION(:), INTENT(in   ) ::   phice, paice  ! ice thickness and concentration
333      REAL(wp), DIMENSION(:), INTENT(inout) ::   pg0, pg1      ! coefficients in linear equation for g(eta)
334      REAL(wp), DIMENSION(:), INTENT(inout) ::   phL, phR      ! min and max value of range over which g(h) > 0
335      !
336      INTEGER  ::   ji           ! horizontal indices
337      REAL(wp) ::   z1_3 , z2_3  ! 1/3 , 2/3
338      REAL(wp) ::   zh13         ! HbL + 1/3 * (HbR - HbL)
339      REAL(wp) ::   zh23         ! HbL + 2/3 * (HbR - HbL)
340      REAL(wp) ::   zdhr         ! 1 / (hR - hL)
341      REAL(wp) ::   zwk1, zwk2   ! temporary variables
342      !!------------------------------------------------------------------
343      !
344      z1_3 = 1._wp / 3._wp
345      z2_3 = 2._wp / 3._wp
346      !
347      DO ji = 1, nidx
348         !
349         IF( paice(ji) > epsi10  .AND. phice(ji) > 0._wp )  THEN
350            !
351            ! Initialize hL and hR
352            phL(ji) = HbL(ji)
353            phR(ji) = HbR(ji)
354            !
355            ! Change hL or hR if hice falls outside central third of range,
356            ! so that hice is in the central third of the range [HL HR]
357            zh13 = z1_3 * ( 2._wp * phL(ji) +         phR(ji) )
358            zh23 = z1_3 * (         phL(ji) + 2._wp * phR(ji) )
359            !
360            IF    ( phice(ji) < zh13 ) THEN   ;   phR(ji) = 3._wp * phice(ji) - 2._wp * phL(ji) ! move HR to the left
361            ELSEIF( phice(ji) > zh23 ) THEN   ;   phL(ji) = 3._wp * phice(ji) - 2._wp * phR(ji) ! move HL to the right
362            ENDIF
363            !
364            ! Compute coefficients of g(eta) = g0 + g1*eta
365            zdhr = 1._wp / (phR(ji) - phL(ji))
366            zwk1 = 6._wp * paice(ji) * zdhr
367            zwk2 = ( phice(ji) - phL(ji) ) * zdhr
368            pg0(ji) = zwk1 * ( z2_3 - zwk2 )                    ! Eq. 14
369            pg1(ji) = 2._wp * zdhr * zwk1 * ( zwk2 - 0.5_wp )   ! Eq. 14
370            !
371         ELSE  ! remap_flag = .false. or a_i < epsi10
372            phL(ji) = 0._wp
373            phR(ji) = 0._wp
374            pg0(ji) = 0._wp
375            pg1(ji) = 0._wp
376         ENDIF
377         !
378      END DO
379      !
380   END SUBROUTINE ice_itd_glinear
381
382
383   SUBROUTINE ice_itd_shiftice( kdonor, pdaice, pdvice )
384      !!------------------------------------------------------------------
385      !!                ***  ROUTINE ice_itd_shiftice ***
386      !!
387      !! ** Purpose :   shift ice across category boundaries, conserving everything
388      !!              ( area, volume, energy, age*vol, and mass of salt )
389      !!------------------------------------------------------------------
390      INTEGER , DIMENSION(:,:), INTENT(in) ::   kdonor   ! donor category index
391      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pdaice   ! ice area transferred across boundary
392      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pdvice   ! ice volume transferred across boundary
393      !
394      INTEGER  ::   ji, jj, jl, jk     ! dummy loop indices
395      INTEGER  ::   ii, ij, jl2, jl1   ! local integers
396      REAL(wp) ::   ztrans             ! ice/snow transferred
397      REAL(wp), DIMENSION(jpij)     ::   zworka, zworkv   ! workspace
398      REAL(wp), DIMENSION(jpij,jpl) ::   zaTsfn           !  -    -
399      !!------------------------------------------------------------------
400         
401      CALL tab_3d_2d( nidx, idxice(1:nidx), ht_i_2d (1:nidx,1:jpl), ht_i   )
402      CALL tab_3d_2d( nidx, idxice(1:nidx), a_i_2d  (1:nidx,1:jpl), a_i    )
403      CALL tab_3d_2d( nidx, idxice(1:nidx), v_i_2d  (1:nidx,1:jpl), v_i    )
404      CALL tab_3d_2d( nidx, idxice(1:nidx), v_s_2d  (1:nidx,1:jpl), v_s    )
405      CALL tab_3d_2d( nidx, idxice(1:nidx), oa_i_2d (1:nidx,1:jpl), oa_i   )
406      CALL tab_3d_2d( nidx, idxice(1:nidx), smv_i_2d(1:nidx,1:jpl), smv_i  )
407      CALL tab_3d_2d( nidx, idxice(1:nidx), a_ip_2d (1:nidx,1:jpl), a_ip   )
408      CALL tab_3d_2d( nidx, idxice(1:nidx), v_ip_2d (1:nidx,1:jpl), v_ip   )
409      CALL tab_3d_2d( nidx, idxice(1:nidx), t_su_2d (1:nidx,1:jpl), t_su   )
410
411      !----------------------------------------------------------------------------------------------
412      ! 1) Define a variable equal to a_i*T_su
413      !----------------------------------------------------------------------------------------------
414      DO jl = 1, jpl
415         DO ji = 1, nidx
416            zaTsfn(ji,jl) = a_i_2d(ji,jl) * t_su_2d(ji,jl)
417         END DO
418      END DO
419     
420      !-------------------------------------------------------------------------------
421      ! 2) Transfer volume and energy between categories
422      !-------------------------------------------------------------------------------
423      DO jl = 1, jpl - 1
424         DO ji = 1, nidx
425            !
426            jl1 = kdonor(ji,jl)
427            !
428            IF( jl1 > 0 ) THEN
429               !
430               IF ( jl1 == jl  ) THEN   ;   jl2 = jl1+1
431               ELSE                     ;   jl2 = jl 
432               ENDIF
433               !
434               IF( v_i_2d(ji,jl1) >= epsi10 ) THEN   ;   zworkv(ji) = pdvice(ji,jl) / v_i_2d(ji,jl1)
435               ELSE                                  ;   zworkv(ji) = 0._wp
436               ENDIF
437               IF( a_i_2d(ji,jl1) >= epsi10 ) THEN   ;   zworka(ji) = pdaice(ji,jl) / a_i_2d(ji,jl1)
438               ELSE                                  ;   zworka(ji) = 0._wp
439               ENDIF
440               !
441               a_i_2d(ji,jl1) = a_i_2d(ji,jl1) - pdaice(ji,jl)       ! Ice areas
442               a_i_2d(ji,jl2) = a_i_2d(ji,jl2) + pdaice(ji,jl)
443               !
444               v_i_2d(ji,jl1) = v_i_2d(ji,jl1) - pdvice(ji,jl)       ! Ice volumes
445               v_i_2d(ji,jl2) = v_i_2d(ji,jl2) + pdvice(ji,jl)
446               !
447               ztrans         = v_s_2d(ji,jl1) * zworkv(ji)          ! Snow volumes
448               v_s_2d(ji,jl1) = v_s_2d(ji,jl1) - ztrans
449               v_s_2d(ji,jl2) = v_s_2d(ji,jl2) + ztrans 
450               !         
451               !                                                     ! Ice age
452               ztrans             = oa_i_2d(ji,jl1) * pdaice(ji,jl)  !!clem: should be * zworka(ji) but it does not work ????
453               oa_i_2d(ji,jl1)    = oa_i_2d(ji,jl1) - ztrans
454               oa_i_2d(ji,jl2)    = oa_i_2d(ji,jl2) + ztrans
455               !
456               ztrans             = smv_i_2d(ji,jl1) * zworkv(ji)    ! Ice salinity
457               !
458               smv_i_2d(ji,jl1)   = smv_i_2d(ji,jl1) - ztrans
459               smv_i_2d(ji,jl2)   = smv_i_2d(ji,jl2) + ztrans
460               !
461               !                                                     ! Surface temperature
462               ztrans          = t_su_2d(ji,jl1) * pdaice(ji,jl)     !!clem: should be * zworka(ji) but it does not work ????
463               zaTsfn(ji,jl1)  = zaTsfn(ji,jl1) - ztrans
464               zaTsfn(ji,jl2)  = zaTsfn(ji,jl2) + ztrans
465               
466               ! MV MP 2016
467               IF ( nn_pnd_scheme > 0 ) THEN
468                  !                                                  ! Pond fraction
469                  ztrans          = a_ip_2d(ji,jl1) * pdaice(ji,jl) !!clem: should be * zworka(ji) but it does not work
470                  a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - ztrans
471                  a_ip_2d(ji,jl2) = a_ip_2d(ji,jl2) + ztrans
472                  !                                                  ! Pond volume (also proportional to da/a)
473                  ztrans          = v_ip_2d(ji,jl1) * pdaice(ji,jl) !!clem: should be * zworka(ji) but it does not work
474                  v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - ztrans
475                  v_ip_2d(ji,jl2) = v_ip_2d(ji,jl2) + ztrans
476               ENDIF
477               ! END MV MP 2016
478               !
479            ENDIF   ! jl1 >0
480         END DO
481         !
482         DO jk = 1, nlay_s         !--- Snow heat content
483            !
484            DO ji = 1, nidx
485               ii = MOD( idxice(ji) - 1, jpi ) + 1
486               ij = ( idxice(ji) - 1 ) / jpi + 1
487               !
488               jl1 = kdonor(ji,jl)
489               !
490               IF( jl1 > 0 ) THEN
491                  IF(jl1 == jl) THEN  ;  jl2 = jl+1
492                  ELSE                ;  jl2 = jl
493                  ENDIF
494                  !
495                  ztrans            = e_s(ii,ij,jk,jl1) * zworkv(ji)
496                  e_s(ii,ij,jk,jl1) = e_s(ii,ij,jk,jl1) - ztrans
497                  e_s(ii,ij,jk,jl2) = e_s(ii,ij,jk,jl2) + ztrans
498               ENDIF
499            END DO
500         END DO
501
502         DO jk = 1, nlay_i         !--- Ice heat content
503            DO ji = 1, nidx
504               ii = MOD( idxice(ji) - 1, jpi ) + 1
505               ij = ( idxice(ji) - 1 ) / jpi + 1
506               !
507               jl1 = kdonor(ji,jl)
508               !
509               IF( jl1 > 0 ) THEN
510                  IF(jl1 == jl) THEN  ;  jl2 = jl+1
511                  ELSE                ;  jl2 = jl
512                  ENDIF
513                  !
514                  ztrans            = e_i(ii,ij,jk,jl1) * zworkv(ji)
515                  e_i(ii,ij,jk,jl1) = e_i(ii,ij,jk,jl1) - ztrans
516                  e_i(ii,ij,jk,jl2) = e_i(ii,ij,jk,jl2) + ztrans
517               ENDIF
518            END DO
519         END DO
520         !
521      END DO                   ! boundaries, 1 to jpl-1
522     
523      !-------------------------------------------------------------------------------
524      ! 3) Update ice thickness and temperature
525      !-------------------------------------------------------------------------------
526      WHERE( a_i_2d(1:nidx,:) >= epsi20 )
527         ht_i_2d(1:nidx,:)  =  v_i_2d(1:nidx,:) / a_i_2d(1:nidx,:) 
528         t_su_2d(1:nidx,:)  =  zaTsfn(1:nidx,:) / a_i_2d(1:nidx,:) 
529      ELSEWHERE
530         ht_i_2d(1:nidx,:)  = 0._wp
531         t_su_2d(1:nidx,:)  = rt0
532      END WHERE
533      !
534      CALL tab_2d_3d( nidx, idxice(1:nidx), ht_i_2d (1:nidx,1:jpl), ht_i  )
535      CALL tab_2d_3d( nidx, idxice(1:nidx), a_i_2d  (1:nidx,1:jpl), a_i   )
536      CALL tab_2d_3d( nidx, idxice(1:nidx), v_i_2d  (1:nidx,1:jpl), v_i   )
537      CALL tab_2d_3d( nidx, idxice(1:nidx), v_s_2d  (1:nidx,1:jpl), v_s   )
538      CALL tab_2d_3d( nidx, idxice(1:nidx), oa_i_2d (1:nidx,1:jpl), oa_i  )
539      CALL tab_2d_3d( nidx, idxice(1:nidx), smv_i_2d(1:nidx,1:jpl), smv_i )
540      CALL tab_2d_3d( nidx, idxice(1:nidx), a_ip_2d (1:nidx,1:jpl), a_ip  )
541      CALL tab_2d_3d( nidx, idxice(1:nidx), v_ip_2d (1:nidx,1:jpl), v_ip  )
542      CALL tab_2d_3d( nidx, idxice(1:nidx), t_su_2d (1:nidx,1:jpl), t_su  )
543      !
544   END SUBROUTINE ice_itd_shiftice
545   
546
547   SUBROUTINE ice_itd_reb
548      !!------------------------------------------------------------------
549      !!                ***  ROUTINE ice_itd_reb ***
550      !!
551      !! ** Purpose : rebin - rebins thicknesses into defined categories
552      !!
553      !! ** Method  : If a category thickness is out of bounds, shift part (for down to top)
554      !!              or entire (for top to down) area, volume, and energy
555      !!              to the neighboring category
556      !!------------------------------------------------------------------
557      INTEGER ::   ji, jj, jl   ! dummy loop indices
558      !
559      INTEGER , DIMENSION(jpij,jpl-1) ::   jdonor           ! donor category index
560      REAL(wp), DIMENSION(jpij,jpl-1) ::   zdaice, zdvice   ! ice area and volume transferred
561      !!------------------------------------------------------------------
562      !
563      jdonor(:,:) = 0
564      zdaice(:,:) = 0._wp
565      zdvice(:,:) = 0._wp
566      !
567      !                       !---------------------------------------
568      DO jl = 1, jpl-1        ! identify thicknesses that are too big
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            END DO
578         END DO
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         END DO
597         !
598         IF( nidx > 0 ) THEN
599            CALL ice_itd_shiftice( jdonor(1:nidx,:), zdaice(1:nidx,:), zdvice(1:nidx,:) )  ! Shift jl=>jl+1
600            ! Reset shift parameters
601            jdonor(1:nidx,jl) = 0
602            zdaice(1:nidx,jl) = 0._wp
603            zdvice(1:nidx,jl) = 0._wp
604         ENDIF
605         !
606      END DO
607
608      !                       !-----------------------------------------
609      DO jl = jpl-1, 1, -1    ! Identify thicknesses that are too small
610         !                    !-----------------------------------------
611         nidx = 0 ; idxice(:) = 0
612         DO jj = 1, jpj
613            DO ji = 1, jpi
614               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
615                  nidx = nidx + 1
616                  idxice( nidx ) = (jj - 1) * jpi + ji                 
617               ENDIF
618            END DO
619         END DO
620         !
621         CALL tab_2d_1d( nidx, idxice(1:nidx), a_i_1d  (1:nidx), a_i(:,:,jl+1)    ) ! jl+1 is ok
622         CALL tab_2d_1d( nidx, idxice(1:nidx), v_i_1d  (1:nidx), v_i(:,:,jl+1)    ) ! jl+1 is ok
623         DO ji = 1, nidx
624            jdonor(ji,jl) = jl + 1
625            zdaice(ji,jl) = a_i_1d(ji) 
626            zdvice(ji,jl) = v_i_1d(ji)
627         END DO
628         !
629         IF( nidx > 0 ) THEN
630            CALL ice_itd_shiftice( jdonor(1:nidx,:), zdaice(1:nidx,:), zdvice(1:nidx,:) )  ! Shift jl+1=>jl
631            ! Reset shift parameters
632            jdonor(1:nidx,jl) = 0
633            zdaice(1:nidx,jl) = 0._wp
634            zdvice(1:nidx,jl) = 0._wp
635         ENDIF
636         !
637      END DO
638      !
639   END SUBROUTINE ice_itd_reb
640
641   SUBROUTINE ice_itd_init
642      !!------------------------------------------------------------------
643      !!                ***  ROUTINE ice_itd_init ***
644      !!
645      !! ** Purpose :   Initializes the ice thickness distribution
646      !! ** Method  :   ...
647      !! ** input   :   Namelist namiceitd
648      !!-------------------------------------------------------------------
649      INTEGER  ::   jl    ! dummy loop index
650      INTEGER  ::   ios   ! Local integer output status for namelist read
651      REAL(wp) ::   zhmax, znum, zden, zalpha   !   -      -
652      !!
653      NAMELIST/namiceitd/ rn_himean
654      !!------------------------------------------------------------------
655      !
656      REWIND( numnam_ice_ref )      ! Namelist namiceitd in reference namelist : Parameters for ice
657      READ  ( numnam_ice_ref, namiceitd, IOSTAT = ios, ERR = 901)
658901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitd in reference namelist', lwp )
659
660      REWIND( numnam_ice_cfg )      ! Namelist namiceitd in configuration namelist : Parameters for ice
661      READ  ( numnam_ice_cfg, namiceitd, IOSTAT = ios, ERR = 902 )
662902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitd in configuration namelist', lwp )
663      IF(lwm) WRITE ( numoni, namiceitd )
664      !
665      IF(lwp) THEN                  ! control print
666         WRITE(numout,*)
667         WRITE(numout,*) 'ice_itd_init : Initialization of ice cat distribution '
668         WRITE(numout,*) '~~~~~~~~~~~~'
669         WRITE(numout,*) '   Namelist namicerun : '
670         WRITE(numout,*) '      mean ice thickness in the domain               rn_himean = ', rn_himean
671      ENDIF
672      !
673      !-----------------------------------!
674      !  Thickness categories boundaries  !
675      !-----------------------------------!
676      !
677      zalpha = 0.05_wp              ! max of each category (from h^(-alpha) function)
678      zhmax  = 3._wp * rn_himean
679      DO jl = 1, jpl
680         znum = jpl * ( zhmax+1 )**zalpha
681         zden = REAL( jpl-jl , wp ) * ( zhmax + 1._wp )**zalpha + REAL( jl , wp )
682         hi_max(jl) = ( znum / zden )**(1./zalpha) - 1
683      END DO
684      !
685      DO jl = 1, jpl                ! mean thickness by category
686         hi_mean(jl) = ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp
687      END DO
688      !
689      hi_max(jpl) = 99._wp          ! set to a big value to ensure that all ice is thinner than hi_max(jpl)
690      !
691      IF(lwp) WRITE(numout,*)
692      IF(lwp) WRITE(numout,*) '   ===>>>   resulting thickness category boundaries :'
693      IF(lwp) WRITE(numout,*) '            hi_max(:)= ', hi_max(0:jpl)
694      !
695   END SUBROUTINE ice_itd_init
696
697#else
698   !!----------------------------------------------------------------------
699   !!   Default option :         Empty module          NO LIM sea-ice model
700   !!----------------------------------------------------------------------
701#endif
702
703   !!======================================================================
704END MODULE iceitd
Note: See TracBrowser for help on using the repository browser.