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

source: branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90 @ 8586

Last change on this file since 8586 was 8586, checked in by gm, 6 years ago

#1911 (ENHANCE-09): PART I.3 - phasing with branch dev_r8183_ICEMODEL revision 8575

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