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_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90 @ 9012

Last change on this file since 9012 was 8967, checked in by clem, 6 years ago

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