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

source: branches/UKMO/dev_merge_2017_GC_couple_pkg/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90 @ 9677

Last change on this file since 9677 was 9677, checked in by dancopsey, 6 years ago

Strip out SVN keywords.

File size: 35.5 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$
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               IF( a_i_2d(ji,jl) > epsi10 )   zdhice(ji,jl) = h_i_2d(ji,jl) - h_ib_2d(ji,jl)   ! clem: if statement is necessary for repro
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) * zworka(ji)
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               sv_i_2d(ji,jl1) = sv_i_2d(ji,jl1) - ztrans
457               sv_i_2d(ji,jl2) = sv_i_2d(ji,jl2) + ztrans
458               !
459               !                                                     ! Surface temperature
460               ztrans          = zaTsfn(ji,jl1) * zworka(ji)
461               zaTsfn(ji,jl1)  = zaTsfn(ji,jl1) - ztrans
462               zaTsfn(ji,jl2)  = zaTsfn(ji,jl2) + ztrans
463               
464               IF ( ln_pnd_H12 ) THEN
465                  !                                                  ! Pond fraction
466                  ztrans          = a_ip_2d(ji,jl1) * zworka(ji)
467                  a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - ztrans
468                  a_ip_2d(ji,jl2) = a_ip_2d(ji,jl2) + ztrans
469                  !                                                  ! Pond volume (also proportional to da/a)
470                  ztrans          = v_ip_2d(ji,jl1) * zworka(ji)
471                  v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - ztrans
472                  v_ip_2d(ji,jl2) = v_ip_2d(ji,jl2) + ztrans
473               ENDIF
474               !
475            ENDIF   ! jl1 >0
476         END DO
477         !
478         DO jk = 1, nlay_s         !--- Snow heat content
479            !
480            DO ji = 1, npti
481               ii = MOD( nptidx(ji) - 1, jpi ) + 1
482               ij = ( nptidx(ji) - 1 ) / jpi + 1
483               !
484               jl1 = kdonor(ji,jl)
485               !
486               IF( jl1 > 0 ) THEN
487                  IF(jl1 == jl) THEN  ;  jl2 = jl+1
488                  ELSE                ;  jl2 = jl
489                  ENDIF
490                  !
491                  ztrans            = e_s(ii,ij,jk,jl1) * zworkv(ji)
492                  e_s(ii,ij,jk,jl1) = e_s(ii,ij,jk,jl1) - ztrans
493                  e_s(ii,ij,jk,jl2) = e_s(ii,ij,jk,jl2) + ztrans
494               ENDIF
495            END DO
496         END DO
497         !
498         DO jk = 1, nlay_i         !--- Ice heat content
499            DO ji = 1, npti
500               ii = MOD( nptidx(ji) - 1, jpi ) + 1
501               ij = ( nptidx(ji) - 1 ) / jpi + 1
502               !
503               jl1 = kdonor(ji,jl)
504               !
505               IF( jl1 > 0 ) THEN
506                  IF(jl1 == jl) THEN  ;  jl2 = jl+1
507                  ELSE                ;  jl2 = jl
508                  ENDIF
509                  !
510                  ztrans            = e_i(ii,ij,jk,jl1) * zworkv(ji)
511                  e_i(ii,ij,jk,jl1) = e_i(ii,ij,jk,jl1) - ztrans
512                  e_i(ii,ij,jk,jl2) = e_i(ii,ij,jk,jl2) + ztrans
513               ENDIF
514            END DO
515         END DO
516         !
517      END DO                   ! boundaries, 1 to jpl-1
518     
519      !-------------------------------------------------------------------------------
520      ! 3) Update ice thickness and temperature
521      !-------------------------------------------------------------------------------
522      WHERE( a_i_2d(1:npti,:) >= epsi20 )
523         h_i_2d(1:npti,:)  =  v_i_2d(1:npti,:) / a_i_2d(1:npti,:) 
524         t_su_2d(1:npti,:)  =  zaTsfn(1:npti,:) / a_i_2d(1:npti,:) 
525      ELSEWHERE
526         h_i_2d(1:npti,:)  = 0._wp
527         t_su_2d(1:npti,:)  = rt0
528      END WHERE
529      !
530      CALL tab_2d_3d( npti, nptidx(1:npti), h_i_2d (1:npti,1:jpl), h_i  )
531      CALL tab_2d_3d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i  )
532      CALL tab_2d_3d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i  )
533      CALL tab_2d_3d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s  )
534      CALL tab_2d_3d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i )
535      CALL tab_2d_3d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i )
536      CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip )
537      CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip )
538      CALL tab_2d_3d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su )
539      !
540   END SUBROUTINE itd_shiftice
541   
542
543   SUBROUTINE ice_itd_reb( kt )
544      !!------------------------------------------------------------------
545      !!                ***  ROUTINE ice_itd_reb ***
546      !!
547      !! ** Purpose : rebin - rebins thicknesses into defined categories
548      !!
549      !! ** Method  : If a category thickness is out of bounds, shift part (for down to top)
550      !!              or entire (for top to down) area, volume, and energy
551      !!              to the neighboring category
552      !!------------------------------------------------------------------
553      INTEGER , INTENT (in) ::   kt      ! Ocean time step
554      INTEGER ::   ji, jj, jl   ! dummy loop indices
555      !
556      INTEGER , DIMENSION(jpij,jpl-1) ::   jdonor           ! donor category index
557      REAL(wp), DIMENSION(jpij,jpl-1) ::   zdaice, zdvice   ! ice area and volume transferred
558      !!------------------------------------------------------------------
559      !
560      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_itd_reb: rebining ice thickness distribution' 
561      !
562      jdonor(:,:) = 0
563      zdaice(:,:) = 0._wp
564      zdvice(:,:) = 0._wp
565      !
566      !                       !---------------------------------------
567      DO jl = 1, jpl-1        ! identify thicknesses that are too big
568         !                    !---------------------------------------
569         npti = 0   ;   nptidx(:) = 0
570         DO jj = 1, jpj
571            DO ji = 1, jpi
572               IF( a_i(ji,jj,jl) > epsi10 .AND. v_i(ji,jj,jl) > (a_i(ji,jj,jl) * hi_max(jl)) ) THEN
573                  npti = npti + 1
574                  nptidx( npti ) = (jj - 1) * jpi + ji                 
575               ENDIF
576            END DO
577         END DO
578         !
579!!clem   CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i(:,:,jl)   )
580         CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d  (1:npti), a_i(:,:,jl)    )
581         CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d  (1:npti), v_i(:,:,jl)    )
582         !
583         DO ji = 1, npti
584            jdonor(ji,jl)  = jl 
585            ! how much of a_i you send in cat sup is somewhat arbitrary
586!!clem: these do not work properly after a restart (I do not know why)
587!!          zdaice(ji,jl)  = a_i_1d(ji) * ( h_i_1d(ji) - hi_max(jl) + epsi10 ) / h_i_1d(ji) 
588!!          zdvice(ji,jl)  = v_i_1d(ji) - ( a_i_1d(ji) - zdaice(ji,jl) ) * ( hi_max(jl) - epsi10 )
589!!clem: these do not work properly after a restart (I do not know why)
590!!            zdaice(ji,jl)  = a_i_1d(ji)
591!!            zdvice(ji,jl)  = v_i_1d(ji)
592!!clem: these are from UCL and work ok
593            zdaice(ji,jl)  = a_i_1d(ji) * 0.5_wp
594            zdvice(ji,jl)  = v_i_1d(ji) - zdaice(ji,jl) * ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp
595         END DO
596         !
597         IF( npti > 0 ) THEN
598            CALL itd_shiftice( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) )  ! Shift jl=>jl+1
599            ! Reset shift parameters
600            jdonor(1:npti,jl) = 0
601            zdaice(1:npti,jl) = 0._wp
602            zdvice(1:npti,jl) = 0._wp
603         ENDIF
604         !
605      END DO
606
607      !                       !-----------------------------------------
608      DO jl = jpl-1, 1, -1    ! Identify thicknesses that are too small
609         !                    !-----------------------------------------
610         npti = 0 ; nptidx(:) = 0
611         DO jj = 1, jpj
612            DO ji = 1, jpi
613               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
614                  npti = npti + 1
615                  nptidx( npti ) = (jj - 1) * jpi + ji                 
616               ENDIF
617            END DO
618         END DO
619         !
620         CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d  (1:npti), a_i(:,:,jl+1)    ) ! jl+1 is ok
621         CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d  (1:npti), v_i(:,:,jl+1)    ) ! jl+1 is ok
622         DO ji = 1, npti
623            jdonor(ji,jl) = jl + 1
624            zdaice(ji,jl) = a_i_1d(ji) 
625            zdvice(ji,jl) = v_i_1d(ji)
626         END DO
627         !
628         IF( npti > 0 ) THEN
629            CALL itd_shiftice( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) )  ! Shift jl+1=>jl
630            ! Reset shift parameters
631            jdonor(1:npti,jl) = 0
632            zdaice(1:npti,jl) = 0._wp
633            zdvice(1:npti,jl) = 0._wp
634         ENDIF
635         !
636      END DO
637      !
638   END SUBROUTINE ice_itd_reb
639
640
641   SUBROUTINE ice_itd_init
642      !!------------------------------------------------------------------
643      !!                ***  ROUTINE ice_itd_init ***
644      !!
645      !! ** Purpose :   Initializes the ice thickness distribution
646      !! ** Method  :   ...
647      !! ** input   :   Namelist namitd
648      !!-------------------------------------------------------------------
649      INTEGER  ::   jl            ! dummy loop index
650      INTEGER  ::   ios, ioptio   ! Local integer output status for namelist read
651      REAL(wp) ::   zhmax, znum, zden, zalpha   !   -      -
652      !
653      NAMELIST/namitd/ ln_cat_hfn, rn_himean, ln_cat_usr, rn_catbnd, rn_himin
654      !!------------------------------------------------------------------
655      !
656      REWIND( numnam_ice_ref )      ! Namelist namitd in reference namelist : Parameters for ice
657      READ  ( numnam_ice_ref, namitd, IOSTAT = ios, ERR = 901)
658901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namitd in reference namelist', lwp )
659      REWIND( numnam_ice_cfg )      ! Namelist namitd in configuration namelist : Parameters for ice
660      READ  ( numnam_ice_cfg, namitd, IOSTAT = ios, ERR = 902 )
661902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namitd in configuration namelist', lwp )
662      IF(lwm) WRITE( numoni, namitd )
663      !
664      IF(lwp) THEN                  ! control print
665         WRITE(numout,*)
666         WRITE(numout,*) 'ice_itd_init: Initialization of ice cat distribution '
667         WRITE(numout,*) '~~~~~~~~~~~~'
668         WRITE(numout,*) '   Namelist namitd: '
669         WRITE(numout,*) '      Ice categories are defined by a function of rn_himean**(-0.05)    ln_cat_hfn = ', ln_cat_hfn
670         WRITE(numout,*) '         mean ice thickness in the domain                               rn_himean  = ', rn_himean
671         WRITE(numout,*) '      Ice categories are defined by rn_catbnd                           ln_cat_usr = ', ln_cat_usr
672         WRITE(numout,*) '      minimum ice thickness                                             rn_himin   = ', rn_himin 
673      ENDIF
674      !
675      !-----------------------------------!
676      !  Thickness categories boundaries  !
677      !-----------------------------------!
678      !                             !== set the choice of ice categories ==!
679      ioptio = 0 
680      IF( ln_cat_hfn ) THEN   ;   ioptio = ioptio + 1   ;   nice_catbnd = np_cathfn    ;   ENDIF
681      IF( ln_cat_usr ) THEN   ;   ioptio = ioptio + 1   ;   nice_catbnd = np_catusr    ;   ENDIF
682      IF( ioptio /= 1 )   CALL ctl_stop( 'ice_itd_init: choose one and only one ice categories boundaries' )
683      !
684      SELECT CASE( nice_catbnd )
685      !                                !------------------------!
686      CASE( np_cathfn )                ! h^(-alpha) function
687         !                             !------------------------!
688         zalpha = 0.05_wp
689         zhmax  = 3._wp * rn_himean
690         hi_max(0) = 0._wp
691         DO jl = 1, jpl
692            znum = jpl * ( zhmax+1 )**zalpha
693            zden = REAL( jpl-jl , wp ) * ( zhmax + 1._wp )**zalpha + REAL( jl , wp )
694            hi_max(jl) = ( znum / zden )**(1./zalpha) - 1
695         END DO
696         !                             !------------------------!
697      CASE( np_catusr )                ! user defined
698         !                             !------------------------!
699         DO jl = 0, jpl
700            hi_max(jl) = rn_catbnd(jl)
701         END DO
702         !
703      END SELECT
704      !
705      DO jl = 1, jpl                ! mean thickness by category
706         hi_mean(jl) = ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp
707      END DO
708      !
709      hi_max(jpl) = 99._wp          ! set to a big value to ensure that all ice is thinner than hi_max(jpl)
710      !
711      IF(lwp) WRITE(numout,*)
712      IF(lwp) WRITE(numout,*) '   ===>>>   resulting thickness category boundaries :'
713      IF(lwp) WRITE(numout,*) '            hi_max(:)= ', hi_max(0:jpl)
714      !
715      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')
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.