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 NEMO/trunk/src/ICE – NEMO

source: NEMO/trunk/src/ICE/iceitd.F90 @ 12808

Last change on this file since 12808 was 12377, checked in by acc, 4 years ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

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