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

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

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

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