source: NEMO/trunk/src/ICE/icedyn_rdgrft.F90 @ 10994

Last change on this file since 10994 was 10994, checked in by clem, 17 months ago

couple of commits to 1) deal with ice concentration reaching 1 (so, no more limitations in the max ice concentration imposed in the namelist), 2) deal with very small negative values that can occur due to roundoff errors

  • Property svn:keywords set to Id
File size: 51.6 KB
Line 
1MODULE icedyn_rdgrft
2   !!======================================================================
3   !!                       ***  MODULE icedyn_rdgrft ***
4   !!    sea-ice : Mechanical impact on ice thickness distribution     
5   !!======================================================================
6   !! History :       !  2006-02  (M. Vancoppenolle) Original code
7   !!            4.0  !  2018     (many people)      SI3 [aka Sea Ice cube]
8   !!----------------------------------------------------------------------
9#if defined key_si3
10   !!----------------------------------------------------------------------
11   !!   'key_si3'                                       SI3 sea-ice model
12   !!----------------------------------------------------------------------
13   !!   ice_dyn_rdgrft       : ridging/rafting of sea ice
14   !!   ice_dyn_rdgrft_init  : initialization of ridging/rafting of sea ice
15   !!   ice_strength         : ice strength calculation
16   !!----------------------------------------------------------------------
17   USE dom_oce        ! ocean domain
18   USE phycst         ! physical constants (ocean directory)
19   USE sbc_oce , ONLY : sss_m, sst_m   ! surface boundary condition: ocean fields
20   USE ice1D          ! sea-ice: thermodynamics
21   USE ice            ! sea-ice: variables
22   USE icetab         ! sea-ice: 1D <==> 2D transformation
23   USE icevar         ! sea-ice: operations
24   USE icectl         ! sea-ice: control prints
25   !
26   USE in_out_manager ! I/O manager
27   USE iom            ! I/O manager library
28   USE lib_mpp        ! MPP library
29   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
30   USE lbclnk         ! lateral boundary conditions (or mpp links)
31   USE timing         ! Timing
32
33   IMPLICIT NONE
34   PRIVATE
35
36   PUBLIC   ice_dyn_rdgrft        ! called by icestp
37   PUBLIC   ice_dyn_rdgrft_init   ! called by icedyn
38   PUBLIC   ice_strength          ! called by icedyn_rhg_evp
39
40   ! Variables shared among ridging subroutines
41   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)     ::   closing_net     ! net rate at which area is removed    (1/s)
42      !                                                               ! (ridging ice area - area of new ridges) / dt
43   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)     ::   opning          ! rate of opening due to divergence/shear
44   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)     ::   closing_gross   ! rate at which area removed, not counting area of new ridges
45   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   apartf          ! participation function; fraction of ridging/closing associated w/ category n
46   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hrmin           ! minimum ridge thickness
47   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hrmax           ! maximum ridge thickness
48   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hraft           ! thickness of rafted ice
49   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hi_hrdg         ! thickness of ridging ice / mean ridge thickness
50   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   aridge          ! participating ice ridging
51   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   araft           ! participating ice rafting
52   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ze_i_2d
53   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ze_s_2d
54   !
55   REAL(wp), PARAMETER ::   hrdg_hi_min = 1.1_wp    ! min ridge thickness multiplier: min(hrdg/hi)
56   REAL(wp), PARAMETER ::   hi_hrft     = 0.5_wp    ! rafting multipliyer: (hi/hraft)
57   !
58   ! ** namelist (namdyn_rdgrft) **
59   LOGICAL  ::   ln_str_H79       ! ice strength parameterization (Hibler79)
60   REAL(wp) ::   rn_pstar         ! determines ice strength, Hibler JPO79
61   REAL(wp) ::   rn_csrdg         ! fraction of shearing energy contributing to ridging           
62   LOGICAL  ::   ln_partf_lin     ! participation function linear (Thorndike et al. (1975))
63   REAL(wp) ::   rn_gstar         !    fractional area of young ice contributing to ridging
64   LOGICAL  ::   ln_partf_exp     ! participation function exponential (Lipscomb et al. (2007))
65   REAL(wp) ::   rn_astar         !    equivalent of G* for an exponential participation function
66   LOGICAL  ::   ln_ridging       ! ridging of ice or not                       
67   REAL(wp) ::   rn_hstar         !    thickness that determines the maximal thickness of ridged ice
68   REAL(wp) ::   rn_porordg       !    initial porosity of ridges (0.3 regular value)
69   REAL(wp) ::   rn_fsnwrdg       !    fractional snow loss to the ocean during ridging
70   REAL(wp) ::   rn_fpndrdg       !    fractional pond loss to the ocean during ridging
71   LOGICAL  ::   ln_rafting       ! rafting of ice or not                       
72   REAL(wp) ::   rn_hraft         !    threshold thickness (m) for rafting / ridging
73   REAL(wp) ::   rn_craft         !    coefficient for smoothness of the hyperbolic tangent in rafting
74   REAL(wp) ::   rn_fsnwrft       !    fractional snow loss to the ocean during rafting
75   REAL(wp) ::   rn_fpndrft       !    fractional pond loss to the ocean during rafting
76   !
77   !!----------------------------------------------------------------------
78   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
79   !! $Id$
80   !! Software governed by the CeCILL licence     (./LICENSE)
81   !!----------------------------------------------------------------------
82CONTAINS
83
84   INTEGER FUNCTION ice_dyn_rdgrft_alloc()
85      !!-------------------------------------------------------------------
86      !!                ***  ROUTINE ice_dyn_rdgrft_alloc ***
87      !!-------------------------------------------------------------------
88      ALLOCATE( closing_net(jpij), opning(jpij)   , closing_gross(jpij),   &
89         &      apartf(jpij,0:jpl), hrmin(jpij,jpl), hraft(jpij,jpl)    , aridge(jpij,jpl),   &
90         &      hrmax(jpij,jpl), hi_hrdg(jpij,jpl)  , araft (jpij,jpl),  &
91         &      ze_i_2d(jpij,nlay_i,jpl), ze_s_2d(jpij,nlay_s,jpl), STAT=ice_dyn_rdgrft_alloc )
92
93      CALL mpp_sum ( 'icedyn_rdgrft', ice_dyn_rdgrft_alloc )
94      IF( ice_dyn_rdgrft_alloc /= 0 )   CALL ctl_stop( 'STOP',  'ice_dyn_rdgrft_alloc: failed to allocate arrays'  )
95      !
96   END FUNCTION ice_dyn_rdgrft_alloc
97
98
99   SUBROUTINE ice_dyn_rdgrft( kt )
100      !!-------------------------------------------------------------------
101      !!                ***  ROUTINE ice_dyn_rdgrft ***
102      !!
103      !! ** Purpose :   computes the mechanical redistribution of ice thickness
104      !!
105      !! ** Method  :   Steps :
106      !!       0) Identify grid cells with ice
107      !!       1) Calculate closing rate, divergence and opening
108      !!       2) Identify grid cells with ridging
109      !!       3) Start ridging iterations
110      !!          - prep = ridged and rafted ice + closing_gross
111      !!          - shift = move ice from one category to another
112      !!
113      !! ** Details
114      !!    step1: The net rate of closing is due to convergence and shear, based on Flato and Hibler (1995).
115      !!           The energy dissipation rate is equal to the net closing rate times the ice strength.
116      !!
117      !!    step3: The gross closing rate is equal to the first two terms (open
118      !!           water closing and thin ice ridging) without the third term
119      !!           (thick, newly ridged ice).
120      !!
121      !! References :   Flato, G. M., and W. D. Hibler III, 1995, JGR, 100, 18,611-18,626.
122      !!                Hibler, W. D. III, 1980, MWR, 108, 1943-1973, 1980.
123      !!                Rothrock, D. A., 1975: JGR, 80, 4514-4519.
124      !!                Thorndike et al., 1975, JGR, 80, 4501-4513.
125      !!                Bitz et al., JGR, 2001
126      !!                Amundrud and Melling, JGR 2005
127      !!                Babko et al., JGR 2002
128      !!
129      !!     This routine is based on CICE code and authors William H. Lipscomb,
130      !!     and Elizabeth C. Hunke, LANL are gratefully acknowledged
131      !!-------------------------------------------------------------------
132      INTEGER, INTENT(in) ::   kt     ! number of iteration
133      !!
134      INTEGER  ::   ji, jj, jk, jl             ! dummy loop index
135      INTEGER  ::   iter, iterate_ridging      ! local integer
136      INTEGER  ::   ipti                       ! local integer
137      REAL(wp) ::   zfac                       ! local scalar
138      INTEGER , DIMENSION(jpij) ::   iptidx        ! compute ridge/raft or not
139      REAL(wp), DIMENSION(jpij) ::   zdivu_adv     ! divu as implied by transport scheme  (1/s)
140      REAL(wp), DIMENSION(jpij) ::   zdivu, zdelt  ! 1D divu_i & delta_i
141      !
142      INTEGER, PARAMETER ::   jp_itermax = 20   
143      !!-------------------------------------------------------------------
144      ! controls
145      IF( ln_timing    )   CALL timing_start('icedyn_rdgrft')                                                             ! timing
146      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icedyn_rdgrft', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
147
148      IF( kt == nit000 ) THEN
149         IF(lwp) WRITE(numout,*)
150         IF(lwp) WRITE(numout,*)'ice_dyn_rdgrft: ice ridging and rafting'
151         IF(lwp) WRITE(numout,*)'~~~~~~~~~~~~~~'
152      ENDIF     
153
154      !--------------------------------
155      ! 0) Identify grid cells with ice
156      !--------------------------------
157      at_i(:,:) = SUM( a_i, dim=3 )
158      !
159      npti = 0   ;   nptidx(:) = 0
160      ipti = 0   ;   iptidx(:) = 0
161      DO jj = 1, jpj
162         DO ji = 1, jpi
163            IF ( at_i(ji,jj) > epsi10 ) THEN
164               npti           = npti + 1
165               nptidx( npti ) = (jj - 1) * jpi + ji
166            ENDIF
167         END DO
168      END DO
169     
170      !--------------------------------------------------------
171      ! 1) Dynamical inputs (closing rate, divergence, opening)
172      !--------------------------------------------------------
173      IF( npti > 0 ) THEN
174       
175         ! just needed here
176         CALL tab_2d_1d( npti, nptidx(1:npti), zdivu   (1:npti)      , divu_i  )
177         CALL tab_2d_1d( npti, nptidx(1:npti), zdelt   (1:npti)      , delta_i )
178         ! needed here and in the iteration loop
179         CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d  (1:npti,1:jpl), a_i   )
180         CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d  (1:npti,1:jpl), v_i   )
181         CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti)      , ato_i )
182
183         DO ji = 1, npti
184            ! closing_net = rate at which open water area is removed + ice area removed by ridging
185            !                                                        - ice area added in new ridges
186            closing_net(ji) = rn_csrdg * 0.5_wp * ( zdelt(ji) - ABS( zdivu(ji) ) ) - MIN( zdivu(ji), 0._wp )
187            !
188            ! divergence given by the advection scheme
189            !   (which may not be equal to divu as computed from the velocity field)
190            IF    ( ln_adv_Pra ) THEN
191               zdivu_adv(ji) = ( 1._wp - ato_i_1d(ji) - SUM( a_i_2d(ji,:) ) ) * r1_rdtice
192            ELSEIF( ln_adv_UMx ) THEN
193               zdivu_adv(ji) = zdivu(ji)
194            ENDIF
195            !
196            IF( zdivu_adv(ji) < 0._wp )   closing_net(ji) = MAX( closing_net(ji), -zdivu_adv(ji) )   ! make sure the closing rate is large enough
197            !                                                                                        ! to give asum = 1.0 after ridging
198            ! Opening rate (non-negative) that will give asum = 1.0 after ridging.
199            opning(ji) = closing_net(ji) + zdivu_adv(ji)
200         END DO
201         !
202         !------------------------------------
203         ! 2) Identify grid cells with ridging
204         !------------------------------------
205         CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d, closing_net )
206
207         DO ji = 1, npti
208            IF( SUM( apartf(ji,1:jpl) ) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
209               ipti = ipti + 1
210               iptidx     (ipti)   = nptidx     (ji)
211               ! adjust to new indices
212               a_i_2d     (ipti,:) = a_i_2d     (ji,:)
213               v_i_2d     (ipti,:) = v_i_2d     (ji,:)
214               ato_i_1d   (ipti)   = ato_i_1d   (ji)
215               closing_net(ipti)   = closing_net(ji)
216               zdivu_adv  (ipti)   = zdivu_adv  (ji)
217               opning     (ipti)   = opning     (ji)
218            ENDIF
219         END DO
220
221      ENDIF
222
223      ! grid cells with ridging
224      nptidx(:) = iptidx(:)
225      npti      = ipti
226
227      !-----------------
228      ! 3) Start ridging
229      !-----------------
230      IF( npti > 0 ) THEN
231         
232         CALL ice_dyn_1d2d( 1 )            ! --- Move to 1D arrays --- !
233
234         iter            = 1
235         iterate_ridging = 1     
236         !                                                        !----------------------!
237         DO WHILE( iterate_ridging > 0 .AND. iter < jp_itermax )  !  ridging iterations  !
238            !                                                     !----------------------!
239            ! Calculate participation function (apartf)
240            !       and transfer      function
241            !       and closing_gross (+correction on opening)
242            CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d, closing_net )
243
244            ! Redistribute area, volume, and energy between categories
245            CALL rdgrft_shift
246
247            ! Do we keep on iterating?
248            !-------------------------
249            ! Check whether a_i + ato_i = 0
250            ! If not, because the closing and opening rates were reduced above, ridge again with new rates
251            iterate_ridging = 0
252            DO ji = 1, npti
253               zfac = 1._wp - ( ato_i_1d(ji) + SUM( a_i_2d(ji,:) ) )
254               IF( ABS( zfac ) < epsi10 ) THEN
255                  closing_net(ji) = 0._wp
256                  opning     (ji) = 0._wp
257                  ato_i_1d   (ji) = MAX( 0._wp, 1._wp - SUM( a_i_2d(ji,:) ) )
258               ELSE
259                  iterate_ridging  = 1
260                  zdivu_adv  (ji) = zfac * r1_rdtice
261                  closing_net(ji) = MAX( 0._wp, -zdivu_adv(ji) )
262                  opning     (ji) = MAX( 0._wp,  zdivu_adv(ji) )
263               ENDIF
264            END DO
265            !
266            iter = iter + 1
267            IF( iter  >  jp_itermax )    CALL ctl_stop( 'STOP',  'icedyn_rdgrft: non-converging ridging scheme'  )
268            !
269         END DO
270
271         CALL ice_dyn_1d2d( 2 )            ! --- Move to 2D arrays --- !
272
273      ENDIF
274   
275      CALL ice_var_agg( 1 ) 
276
277      ! controls
278      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icedyn_rdgrft', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
279      IF( ln_ctl       )   CALL ice_prt3D   ('icedyn_rdgrft')                                                             ! prints
280      IF( ln_timing    )   CALL timing_stop ('icedyn_rdgrft')                                                             ! timing
281      !
282   END SUBROUTINE ice_dyn_rdgrft
283
284
285   SUBROUTINE rdgrft_prep( pa_i, pv_i, pato_i, pclosing_net )
286      !!-------------------------------------------------------------------
287      !!                ***  ROUTINE rdgrft_prep ***
288      !!
289      !! ** Purpose :   preparation for ridging calculations
290      !!
291      !! ** Method  :   Compute the thickness distribution of the ice and open water
292      !!                participating in ridging and of the resulting ridges.
293      !!-------------------------------------------------------------------
294      REAL(wp), DIMENSION(:)  , INTENT(in) ::   pato_i, pclosing_net 
295      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pa_i, pv_i 
296      !!
297      INTEGER  ::   ji, jl                     ! dummy loop indices
298      REAL(wp) ::   z1_gstar, z1_astar, zhmean, zfac   ! local scalar
299      REAL(wp), DIMENSION(jpij)        ::   zasum, z1_asum, zaksum   ! sum of a_i+ato_i and reverse
300      REAL(wp), DIMENSION(jpij,jpl)    ::   zhi                      ! ice thickness
301      REAL(wp), DIMENSION(jpij,-1:jpl) ::   zGsum                    ! zGsum(n) = sum of areas in categories 0 to n
302      !--------------------------------------------------------------------
303
304      z1_gstar = 1._wp / rn_gstar
305      z1_astar = 1._wp / rn_astar
306
307      !                       ! Ice thickness needed for rafting
308      WHERE( pa_i(1:npti,:) > epsi20 )   ;   zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:)
309      ELSEWHERE                          ;   zhi(1:npti,:) = 0._wp
310      END WHERE
311
312      ! 1) Participation function (apartf): a(h) = b(h).g(h)
313      !-----------------------------------------------------------------
314      ! Compute the participation function = total area lost due to ridging/closing
315      ! This is analogous to
316      !   a(h) = b(h)g(h) as defined in Thorndike et al. (1975).
317      !   assuming b(h) = (2/Gstar) * (1 - G(h)/Gstar).
318      !
319      ! apartf = integrating b(h)g(h) between the category boundaries
320      ! apartf is always >= 0 and SUM(apartf(0:jpl))=1
321      !-----------------------------------------------------------------
322      !
323      ! Compute total area of ice plus open water.
324      ! This is in general not equal to one because of divergence during transport
325      zasum(1:npti) = pato_i(1:npti) + SUM( pa_i(1:npti,:), dim=2 )
326      !
327      WHERE( zasum(1:npti) > epsi20 )   ;   z1_asum(1:npti) = 1._wp / zasum(1:npti)
328      ELSEWHERE                         ;   z1_asum(1:npti) = 0._wp
329      END WHERE
330      !
331      ! Compute cumulative thickness distribution function
332      ! Compute the cumulative thickness distribution function zGsum,
333      ! where zGsum(n) is the fractional area in categories 0 to n.
334      ! initial value (in h = 0) = open water area
335      zGsum(1:npti,-1) = 0._wp
336      zGsum(1:npti,0 ) = pato_i(1:npti) * z1_asum(1:npti)
337      DO jl = 1, jpl
338         zGsum(1:npti,jl) = ( pato_i(1:npti) + SUM( pa_i(1:npti,1:jl), dim=2 ) ) * z1_asum(1:npti)  ! sum(1:jl) is ok (and not jpl)
339      END DO
340      !
341      IF( ln_partf_lin ) THEN          !--- Linear formulation (Thorndike et al., 1975)
342         DO jl = 0, jpl   
343            DO ji = 1, npti
344               IF    ( zGsum(ji,jl)   < rn_gstar ) THEN
345                  apartf(ji,jl) = z1_gstar * ( zGsum(ji,jl) - zGsum(ji,jl-1) ) * &
346                     &                       ( 2._wp - ( zGsum(ji,jl-1) + zGsum(ji,jl) ) * z1_gstar )
347               ELSEIF( zGsum(ji,jl-1) < rn_gstar ) THEN
348                  apartf(ji,jl) = z1_gstar * ( rn_gstar     - zGsum(ji,jl-1) ) *  &
349                     &                       ( 2._wp - ( zGsum(ji,jl-1) + rn_gstar        ) * z1_gstar )
350               ELSE
351                  apartf(ji,jl) = 0._wp
352               ENDIF
353            END DO
354         END DO
355         !
356      ELSEIF( ln_partf_exp ) THEN      !--- Exponential, more stable formulation (Lipscomb et al, 2007)
357         !                       
358         zfac = 1._wp / ( 1._wp - EXP(-z1_astar) )
359         DO jl = -1, jpl
360            DO ji = 1, npti
361               zGsum(ji,jl) = EXP( -zGsum(ji,jl) * z1_astar ) * zfac
362            END DO
363         END DO
364         DO jl = 0, jpl
365            DO ji = 1, npti
366               apartf(ji,jl) = zGsum(ji,jl-1) - zGsum(ji,jl)
367            END DO
368         END DO
369         !
370      ENDIF
371
372      !                                !--- Ridging and rafting participation concentrations
373      IF( ln_rafting .AND. ln_ridging ) THEN             !- ridging & rafting
374         DO jl = 1, jpl
375            DO ji = 1, npti
376               aridge(ji,jl) = ( 1._wp + TANH ( rn_craft * ( zhi(ji,jl) - rn_hraft ) ) ) * 0.5_wp * apartf(ji,jl)
377               araft (ji,jl) = apartf(ji,jl) - aridge(ji,jl)
378            END DO
379         END DO
380      ELSEIF( ln_ridging .AND. .NOT. ln_rafting ) THEN   !- ridging alone
381         DO jl = 1, jpl
382            DO ji = 1, npti
383               aridge(ji,jl) = apartf(ji,jl)
384               araft (ji,jl) = 0._wp
385            END DO
386         END DO
387      ELSEIF( ln_rafting .AND. .NOT. ln_ridging ) THEN   !- rafting alone   
388         DO jl = 1, jpl
389            DO ji = 1, npti
390               aridge(ji,jl) = 0._wp
391               araft (ji,jl) = apartf(ji,jl)
392            END DO
393         END DO
394      ELSE                                               !- no ridging & no rafting
395         DO jl = 1, jpl
396            DO ji = 1, npti
397               aridge(ji,jl) = 0._wp
398               araft (ji,jl) = 0._wp         
399            END DO
400         END DO
401      ENDIF
402
403      ! 2) Transfer function
404      !-----------------------------------------------------------------
405      ! Compute max and min ridged ice thickness for each ridging category.
406      ! Assume ridged ice is uniformly distributed between hrmin and hrmax.
407      !
408      ! This parameterization is a modified version of Hibler (1980).
409      ! The mean ridging thickness, zhmean, is proportional to hi^(0.5)
410      !  and for very thick ridging ice must be >= hrdg_hi_min*hi
411      !
412      ! The minimum ridging thickness, hrmin, is equal to 2*hi
413      !  (i.e., rafting) and for very thick ridging ice is
414      !  constrained by hrmin <= (zhmean + hi)/2.
415      !
416      ! The maximum ridging thickness, hrmax, is determined by zhmean and hrmin.
417      !
418      ! These modifications have the effect of reducing the ice strength
419      ! (relative to the Hibler formulation) when very thick ice is ridging.
420      !
421      ! zaksum = net area removed/ total area removed
422      ! where total area removed = area of ice that ridges
423      !         net area removed = total area removed - area of new ridges
424      !-----------------------------------------------------------------
425      zfac = 1._wp / hi_hrft
426      zaksum(1:npti) = apartf(1:npti,0)
427      !
428      DO jl = 1, jpl
429         DO ji = 1, npti
430            IF ( apartf(ji,jl) > 0._wp ) THEN
431               zhmean         = MAX( SQRT( rn_hstar * zhi(ji,jl) ), zhi(ji,jl) * hrdg_hi_min )
432               hrmin  (ji,jl) = MIN( 2._wp * zhi(ji,jl), 0.5_wp * ( zhmean + zhi(ji,jl) ) )
433               hrmax  (ji,jl) = 2._wp * zhmean - hrmin(ji,jl)
434               hraft  (ji,jl) = zhi(ji,jl) * zfac
435               hi_hrdg(ji,jl) = zhi(ji,jl) / MAX( zhmean, epsi20 )
436               !
437               ! Normalization factor : zaksum, ensures mass conservation
438               zaksum(ji) = zaksum(ji) + aridge(ji,jl) * ( 1._wp - hi_hrdg(ji,jl) )    &
439                  &                    + araft (ji,jl) * ( 1._wp - hi_hrft )
440            ELSE
441               hrmin  (ji,jl) = 0._wp 
442               hrmax  (ji,jl) = 0._wp 
443               hraft  (ji,jl) = 0._wp 
444               hi_hrdg(ji,jl) = 1._wp
445            ENDIF
446         END DO
447      END DO
448      !
449      ! 3) closing_gross
450      !-----------------
451      ! Based on the ITD of ridging and ridged ice, convert the net closing rate to a gross closing rate. 
452      ! NOTE: 0 < aksum <= 1
453      WHERE( zaksum(1:npti) > epsi20 )   ;   closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti)
454      ELSEWHERE                          ;   closing_gross(1:npti) = 0._wp
455      END WHERE
456     
457      ! correction to closing rate if excessive ice removal
458      !----------------------------------------------------
459      ! Reduce the closing rate if more than 100% of any ice category would be removed
460      ! Reduce the opening rate in proportion
461      DO jl = 1, jpl
462         DO ji = 1, npti
463            zfac = apartf(ji,jl) * closing_gross(ji) * rdt_ice
464            IF( zfac > pa_i(ji,jl) .AND. apartf(ji,jl) /= 0._wp ) THEN
465               closing_gross(ji) = pa_i(ji,jl) / apartf(ji,jl) * r1_rdtice
466            ENDIF
467         END DO
468      END DO     
469
470      ! 4) correction to opening if excessive open water removal
471      !---------------------------------------------------------
472      ! Reduce the closing rate if more than 100% of the open water would be removed
473      ! Reduce the opening rate in proportion
474      DO ji = 1, npti 
475         zfac = pato_i(ji) + ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * rdt_ice
476         IF( zfac < 0._wp ) THEN           ! would lead to negative ato_i
477            opning(ji) = apartf(ji,0) * closing_gross(ji) - pato_i(ji) * r1_rdtice 
478         ELSEIF( zfac > zasum(ji) ) THEN   ! would lead to ato_i > asum
479            opning(ji) = apartf(ji,0) * closing_gross(ji) + ( zasum(ji) - pato_i(ji) ) * r1_rdtice 
480         ENDIF
481      END DO
482      !
483   END SUBROUTINE rdgrft_prep
484
485
486   SUBROUTINE rdgrft_shift
487      !!-------------------------------------------------------------------
488      !!                ***  ROUTINE rdgrft_shift ***
489      !!
490      !! ** Purpose :   shift ridging ice among thickness categories of ice thickness
491      !!
492      !! ** Method  :   Remove area, volume, and energy from each ridging category
493      !!                and add to thicker ice categories.
494      !!-------------------------------------------------------------------
495      !
496      INTEGER  ::   ji, jj, jl, jl1, jl2, jk   ! dummy loop indices
497      REAL(wp) ::   hL, hR, farea              ! left and right limits of integration and new area going to jl2
498      REAL(wp) ::   vsw                        ! vol of water trapped into ridges
499      REAL(wp) ::   afrdg, afrft               ! fraction of category area ridged/rafted
500      REAL(wp)                  ::   airdg1, oirdg1, aprdg1, virdg1, sirdg1
501      REAL(wp)                  ::   airft1, oirft1, aprft1
502      REAL(wp), DIMENSION(jpij) ::   airdg2, oirdg2, aprdg2, virdg2, sirdg2, vsrdg, vprdg  ! area etc of new ridges
503      REAL(wp), DIMENSION(jpij) ::   airft2, oirft2, aprft2, virft , sirft , vsrft, vprft  ! area etc of rafted ice
504      !
505      REAL(wp), DIMENSION(jpij) ::   ersw             ! enth of water trapped into ridges
506      REAL(wp), DIMENSION(jpij) ::   zswitch, fvol    ! new ridge volume going to jl2
507      REAL(wp), DIMENSION(jpij) ::   z1_ai            ! 1 / a
508      REAL(wp), DIMENSION(jpij) ::   zvti             ! sum(v_i)
509      !
510      REAL(wp), DIMENSION(jpij,nlay_s) ::   esrft     ! snow energy of rafting ice
511      REAL(wp), DIMENSION(jpij,nlay_i) ::   eirft     ! ice  energy of rafting ice
512      REAL(wp), DIMENSION(jpij,nlay_s) ::   esrdg     ! enth*volume of new ridges     
513      REAL(wp), DIMENSION(jpij,nlay_i) ::   eirdg     ! enth*volume of new ridges
514      !
515      INTEGER , DIMENSION(jpij) ::   itest_rdg, itest_rft   ! test for conservation
516      !!-------------------------------------------------------------------
517      !
518      zvti(1:npti) = SUM( v_i_2d(1:npti,:), dim=2 )   ! total ice volume
519      !
520      ! 1) Change in open water area due to closing and opening
521      !--------------------------------------------------------
522      DO ji = 1, npti
523         ato_i_1d(ji) = MAX( 0._wp, ato_i_1d(ji) + ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * rdt_ice )
524      END DO
525     
526      ! 2) compute categories in which ice is removed (jl1)
527      !----------------------------------------------------
528      DO jl1 = 1, jpl
529
530         CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d(1:npti), s_i(:,:,jl1) )
531
532         DO ji = 1, npti
533
534            IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN   ! only if ice is ridging
535
536               IF( a_i_2d(ji,jl1) > epsi20 ) THEN   ;   z1_ai(ji) = 1._wp / a_i_2d(ji,jl1)
537               ELSE                                 ;   z1_ai(ji) = 0._wp
538               ENDIF
539               
540               ! area of ridging / rafting ice (airdg1) and of new ridge (airdg2)
541               airdg1 = aridge(ji,jl1) * closing_gross(ji) * rdt_ice
542               airft1 = araft (ji,jl1) * closing_gross(ji) * rdt_ice
543
544               airdg2(ji) = airdg1 * hi_hrdg(ji,jl1)
545               airft2(ji) = airft1 * hi_hrft
546
547               ! ridging /rafting fractions
548               afrdg = airdg1 * z1_ai(ji)
549               afrft = airft1 * z1_ai(ji)
550
551               ! volume and enthalpy (J/m2, >0) of seawater trapped into ridges
552               IF    ( zvti(ji) <= 10. ) THEN ; vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg                                           ! v <= 10m then porosity = rn_porordg
553               ELSEIF( zvti(ji) >= 20. ) THEN ; vsw = 0._wp                                                                         ! v >= 20m then porosity = 0
554               ELSE                           ; vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg * MAX( 0._wp, 2._wp - 0.1_wp * zvti(ji) ) ! v > 10m and v < 20m then porosity = linear transition to 0
555               ENDIF
556               ersw(ji) = -rhoi * vsw * rcp * sst_1d(ji)   ! clem: if sst>0, then ersw <0 (is that possible?)
557
558               ! volume etc of ridging / rafting ice and new ridges (vi, vs, sm, oi, es, ei)
559               virdg1     = v_i_2d (ji,jl1)   * afrdg
560               virdg2(ji) = v_i_2d (ji,jl1)   * afrdg + vsw
561               vsrdg(ji)  = v_s_2d (ji,jl1)   * afrdg
562               sirdg1     = sv_i_2d(ji,jl1)   * afrdg
563               sirdg2(ji) = sv_i_2d(ji,jl1)   * afrdg + vsw * sss_1d(ji)
564               oirdg1     = oa_i_2d(ji,jl1)   * afrdg
565               oirdg2(ji) = oa_i_2d(ji,jl1)   * afrdg * hi_hrdg(ji,jl1) 
566
567               virft(ji)  = v_i_2d (ji,jl1)   * afrft
568               vsrft(ji)  = v_s_2d (ji,jl1)   * afrft
569               sirft(ji)  = sv_i_2d(ji,jl1)   * afrft 
570               oirft1     = oa_i_2d(ji,jl1)   * afrft 
571               oirft2(ji) = oa_i_2d(ji,jl1)   * afrft * hi_hrft 
572
573               IF ( ln_pnd_H12 ) THEN
574                  aprdg1     = a_ip_2d(ji,jl1) * afrdg
575                  aprdg2(ji) = a_ip_2d(ji,jl1) * afrdg * hi_hrdg(ji,jl1)
576                  vprdg (ji) = v_ip_2d(ji,jl1) * afrdg
577                  aprft1     = a_ip_2d(ji,jl1) * afrft
578                  aprft2(ji) = a_ip_2d(ji,jl1) * afrft * hi_hrft
579                  vprft (ji) = v_ip_2d(ji,jl1) * afrft
580               ENDIF
581
582               ! Ice-ocean exchanges associated with ice porosity
583               wfx_dyn_1d(ji) = wfx_dyn_1d(ji) - vsw * rhoi * r1_rdtice   ! increase in ice volume due to seawater frozen in voids
584               sfx_dyn_1d(ji) = sfx_dyn_1d(ji) - vsw * sss_1d(ji) * rhoi * r1_rdtice
585               hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ersw(ji) * r1_rdtice          ! > 0 [W.m-2]
586
587               ! Put the snow lost by ridging into the ocean
588               !  Note that esrdg > 0; the ocean must cool to melt snow. If the ocean temp = Tf already, new ice must grow.
589               wfx_snw_dyn_1d(ji) = wfx_snw_dyn_1d(ji) + ( rhos * vsrdg(ji) * ( 1._wp - rn_fsnwrdg )   &   ! fresh water source for ocean
590                  &                                      + rhos * vsrft(ji) * ( 1._wp - rn_fsnwrft ) ) * r1_rdtice
591
592               ! virtual salt flux to keep salinity constant
593               IF( nn_icesal /= 2 )  THEN
594                  sirdg2(ji)     = sirdg2(ji)     - vsw * ( sss_1d(ji) - s_i_1d(ji) )        ! ridge salinity = s_i
595                  sfx_bri_1d(ji) = sfx_bri_1d(ji) + sss_1d(ji) * vsw * rhoi * r1_rdtice  &  ! put back sss_m into the ocean
596                     &                            - s_i_1d(ji) * vsw * rhoi * r1_rdtice     ! and get  s_i  from the ocean
597               ENDIF
598
599               ! Remove area, volume of new ridge to each category jl1
600               !------------------------------------------------------
601               a_i_2d (ji,jl1) = a_i_2d (ji,jl1) - airdg1    - airft1
602               v_i_2d (ji,jl1) = v_i_2d (ji,jl1) - virdg1    - virft(ji)
603               v_s_2d (ji,jl1) = v_s_2d (ji,jl1) - vsrdg(ji) - vsrft(ji)
604               sv_i_2d(ji,jl1) = sv_i_2d(ji,jl1) - sirdg1    - sirft(ji)
605               oa_i_2d(ji,jl1) = oa_i_2d(ji,jl1) - oirdg1    - oirft1
606               IF ( ln_pnd_H12 ) THEN
607                  a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - aprdg1    - aprft1
608                  v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - vprdg(ji) - vprft(ji)
609               ENDIF
610            ENDIF
611
612         END DO ! ji
613
614         ! special loop for e_s because of layers jk
615         DO jk = 1, nlay_s
616            DO ji = 1, npti
617               IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
618                  ! Compute ridging /rafting fractions
619                  afrdg = aridge(ji,jl1) * closing_gross(ji) * rdt_ice * z1_ai(ji)
620                  afrft = araft (ji,jl1) * closing_gross(ji) * rdt_ice * z1_ai(ji)
621                  ! Compute ridging /rafting ice and new ridges for es
622                  esrdg(ji,jk) = ze_s_2d (ji,jk,jl1) * afrdg
623                  esrft(ji,jk) = ze_s_2d (ji,jk,jl1) * afrft
624                  ! Put the snow lost by ridging into the ocean
625                  hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ( - esrdg(ji,jk) * ( 1._wp - rn_fsnwrdg )   &                 ! heat sink for ocean (<0, W.m-2)
626                     &                                - esrft(ji,jk) * ( 1._wp - rn_fsnwrft ) ) * r1_rdtice
627                  !
628                  ! Remove energy of new ridge to each category jl1
629                  !-------------------------------------------------
630                  ze_s_2d(ji,jk,jl1) = ze_s_2d(ji,jk,jl1) * ( 1._wp - afrdg - afrft ) 
631               ENDIF
632            END DO
633         END DO
634                 
635         ! special loop for e_i because of layers jk
636         DO jk = 1, nlay_i
637            DO ji = 1, npti
638               IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
639                  ! Compute ridging /rafting fractions
640                  afrdg = aridge(ji,jl1) * closing_gross(ji) * rdt_ice * z1_ai(ji)
641                  afrft = araft (ji,jl1) * closing_gross(ji) * rdt_ice * z1_ai(ji)
642                  ! Compute ridging ice and new ridges for ei
643                  eirdg(ji,jk) = ze_i_2d (ji,jk,jl1) * afrdg + ersw(ji) * r1_nlay_i
644                  eirft(ji,jk) = ze_i_2d (ji,jk,jl1) * afrft
645                  !
646                  ! Remove energy of new ridge to each category jl1
647                  !-------------------------------------------------
648                  ze_i_2d(ji,jk,jl1) = ze_i_2d(ji,jk,jl1) * ( 1._wp - afrdg - afrft ) 
649               ENDIF
650            END DO
651         END DO
652         
653         ! 3) compute categories in which ice is added (jl2)
654         !--------------------------------------------------
655         itest_rdg(1:npti) = 0
656         itest_rft(1:npti) = 0
657         DO jl2  = 1, jpl 
658            !
659            DO ji = 1, npti
660
661               IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
662
663                  ! Compute the fraction of ridged ice area and volume going to thickness category jl2
664                  IF( hrmin(ji,jl1) <= hi_max(jl2) .AND. hrmax(ji,jl1) > hi_max(jl2-1) ) THEN
665                     hL = MAX( hrmin(ji,jl1), hi_max(jl2-1) )
666                     hR = MIN( hrmax(ji,jl1), hi_max(jl2)   )
667                     farea    = ( hR      - hL      ) / ( hrmax(ji,jl1)                 - hrmin(ji,jl1)                 )
668                     fvol(ji) = ( hR * hR - hL * hL ) / ( hrmax(ji,jl1) * hrmax(ji,jl1) - hrmin(ji,jl1) * hrmin(ji,jl1) )
669                     !
670                     itest_rdg(ji) = 1   ! test for conservation
671                  ELSE
672                     farea    = 0._wp 
673                     fvol(ji) = 0._wp                 
674                  ENDIF
675
676                  ! Compute the fraction of rafted ice area and volume going to thickness category jl2
677                  IF( hraft(ji,jl1) <= hi_max(jl2) .AND. hraft(ji,jl1) >  hi_max(jl2-1) ) THEN
678                     zswitch(ji) = 1._wp
679                     !
680                     itest_rft(ji) = 1   ! test for conservation
681                  ELSE
682                     zswitch(ji) = 0._wp
683                  ENDIF
684                  !
685                  ! Patch to ensure perfect conservation if ice thickness goes mad
686                  ! Sometimes thickness is larger than hi_max(jpl) because of advection scheme (for very small areas)
687                  ! Then ice volume is removed from one category but the ridging/rafting scheme
688                  ! does not know where to move it, leading to a conservation issue. 
689                  IF( itest_rdg(ji) == 0 .AND. jl2 == jpl ) THEN   ;   farea = 1._wp   ;   fvol(ji) = 1._wp   ;   ENDIF
690                  IF( itest_rft(ji) == 0 .AND. jl2 == jpl )      zswitch(ji) = 1._wp
691                  !
692                  ! Add area, volume of new ridge to category jl2
693                  !----------------------------------------------
694                  a_i_2d (ji,jl2) = a_i_2d (ji,jl2) + ( airdg2(ji) * farea    + airft2(ji) * zswitch(ji) )
695                  oa_i_2d(ji,jl2) = oa_i_2d(ji,jl2) + ( oirdg2(ji) * farea    + oirft2(ji) * zswitch(ji) )
696                  v_i_2d (ji,jl2) = v_i_2d (ji,jl2) + ( virdg2(ji) * fvol(ji) + virft (ji) * zswitch(ji) )
697                  sv_i_2d(ji,jl2) = sv_i_2d(ji,jl2) + ( sirdg2(ji) * fvol(ji) + sirft (ji) * zswitch(ji) )
698                  v_s_2d (ji,jl2) = v_s_2d (ji,jl2) + ( vsrdg (ji) * rn_fsnwrdg * fvol(ji)  +  &
699                     &                                  vsrft (ji) * rn_fsnwrft * zswitch(ji) )
700                  IF ( ln_pnd_H12 ) THEN
701                     v_ip_2d (ji,jl2) = v_ip_2d(ji,jl2) + (   vprdg (ji) * rn_fpndrdg * fvol   (ji)   &
702                        &                                   + vprft (ji) * rn_fpndrft * zswitch(ji)   )
703                     a_ip_2d (ji,jl2) = a_ip_2d(ji,jl2) + (   aprdg2(ji) * rn_fpndrdg * farea         & 
704                        &                                   + aprft2(ji) * rn_fpndrft * zswitch(ji)   )
705                  ENDIF
706                 
707               ENDIF
708
709            END DO
710            ! Add snow energy of new ridge to category jl2
711            !---------------------------------------------
712            DO jk = 1, nlay_s
713               DO ji = 1, npti
714                  IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp )   &
715                     &   ze_s_2d(ji,jk,jl2) = ze_s_2d(ji,jk,jl2) + ( esrdg(ji,jk) * rn_fsnwrdg * fvol(ji)  +  &
716                     &                                               esrft(ji,jk) * rn_fsnwrft * zswitch(ji) )
717               END DO
718            END DO
719            ! Add ice energy of new ridge to category jl2
720            !--------------------------------------------
721            DO jk = 1, nlay_i
722               DO ji = 1, npti
723                  IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp )   &
724                     &   ze_i_2d(ji,jk,jl2) = ze_i_2d(ji,jk,jl2) + eirdg(ji,jk) * fvol(ji) + eirft(ji,jk) * zswitch(ji)                 
725               END DO
726            END DO
727            !
728         END DO ! jl2
729         !
730      END DO ! jl1
731      !
732      ! roundoff errors
733      !----------------
734      ! In case ridging/rafting lead to very small negative values (sometimes it happens)
735      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 )
736      !
737   END SUBROUTINE rdgrft_shift
738
739
740   SUBROUTINE ice_strength
741      !!----------------------------------------------------------------------
742      !!                ***  ROUTINE ice_strength ***
743      !!
744      !! ** Purpose :   computes ice strength used in dynamics routines of ice thickness
745      !!
746      !! ** Method  :   Compute the strength of the ice pack, defined as the energy (J m-2)
747      !!              dissipated per unit area removed from the ice pack under compression,
748      !!              and assumed proportional to the change in potential energy caused
749      !!              by ridging. Note that only Hibler's formulation is stable and that
750      !!              ice strength has to be smoothed
751      !!----------------------------------------------------------------------
752      INTEGER             ::   ji, jj, jl  ! dummy loop indices
753      INTEGER             ::   ismooth     ! smoothing the resistance to deformation
754      INTEGER             ::   itframe     ! number of time steps for the P smoothing
755      REAL(wp)            ::   zp, z1_3    ! local scalars
756      REAL(wp), DIMENSION(jpi,jpj) ::   zworka           ! temporary array used here
757      REAL(wp), DIMENSION(jpi,jpj) ::   zstrp1, zstrp2   ! strength at previous time steps
758      !!----------------------------------------------------------------------
759      !                              !--------------------------------------------------!
760      IF( ln_str_H79 ) THEN          ! Ice strength => Hibler (1979) method             !
761      !                              !--------------------------------------------------!
762         strength(:,:) = rn_pstar * SUM( v_i(:,:,:), dim=3 ) * EXP( -rn_crhg * ( 1._wp - SUM( a_i(:,:,:), dim=3 ) ) )
763         ismooth = 1
764         !                           !--------------------------------------------------!
765      ELSE                           ! Zero strength                                    !
766         !                           !--------------------------------------------------!
767         strength(:,:) = 0._wp
768         ismooth = 0
769      ENDIF
770      !                              !--------------------------------------------------!
771      SELECT CASE( ismooth )         ! Smoothing ice strength                           !
772      !                              !--------------------------------------------------!
773      CASE( 1 )               !--- Spatial smoothing
774         DO jj = 2, jpjm1
775            DO ji = 2, jpim1
776               IF ( SUM( a_i(ji,jj,:) ) > 0._wp ) THEN
777                  zworka(ji,jj) = ( 4.0 * strength(ji,jj)              &
778                     &                  + strength(ji-1,jj) * tmask(ji-1,jj,1) + strength(ji+1,jj) * tmask(ji+1,jj,1) & 
779                     &                  + strength(ji,jj-1) * tmask(ji,jj-1,1) + strength(ji,jj+1) * tmask(ji,jj+1,1) &
780                     &            ) / ( 4.0 + tmask(ji-1,jj,1) + tmask(ji+1,jj,1) + tmask(ji,jj-1,1) + tmask(ji,jj+1,1) )
781               ELSE
782                  zworka(ji,jj) = 0._wp
783               ENDIF
784            END DO
785         END DO
786         
787         DO jj = 2, jpjm1
788            DO ji = 2, jpim1
789               strength(ji,jj) = zworka(ji,jj)
790            END DO
791         END DO
792         CALL lbc_lnk( 'icedyn_rdgrft', strength, 'T', 1. )
793         !
794      CASE( 2 )               !--- Temporal smoothing
795         IF ( kt_ice == nit000 ) THEN
796            zstrp1(:,:) = 0._wp
797            zstrp2(:,:) = 0._wp
798         ENDIF
799         !
800         DO jj = 2, jpjm1
801            DO ji = 2, jpim1
802               IF ( SUM( a_i(ji,jj,:) ) > 0._wp ) THEN
803                  itframe = 1 ! number of time steps for the running mean
804                  IF ( zstrp1(ji,jj) > 0._wp ) itframe = itframe + 1
805                  IF ( zstrp2(ji,jj) > 0._wp ) itframe = itframe + 1
806                  zp = ( strength(ji,jj) + zstrp1(ji,jj) + zstrp2(ji,jj) ) / itframe
807                  zstrp2  (ji,jj) = zstrp1  (ji,jj)
808                  zstrp1  (ji,jj) = strength(ji,jj)
809                  strength(ji,jj) = zp
810               ENDIF
811            END DO
812         END DO
813         CALL lbc_lnk( 'icedyn_rdgrft', strength, 'T', 1. )
814         !
815      END SELECT
816      !
817   END SUBROUTINE ice_strength
818
819   
820   SUBROUTINE ice_dyn_1d2d( kn )
821      !!-----------------------------------------------------------------------
822      !!                   ***  ROUTINE ice_dyn_1d2d ***
823      !!                 
824      !! ** Purpose :   move arrays from 1d to 2d and the reverse
825      !!-----------------------------------------------------------------------
826      INTEGER, INTENT(in) ::   kn   ! 1= from 2D to 1D   ;   2= from 1D to 2D
827      !
828      INTEGER ::   jl, jk   ! dummy loop indices
829      !!-----------------------------------------------------------------------
830      !
831      SELECT CASE( kn )
832      !                    !---------------------!
833      CASE( 1 )            !==  from 2D to 1D  ==!
834         !                 !---------------------!
835         ! fields used but not modified
836         CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d(1:npti), sss_m(:,:) )
837         CALL tab_2d_1d( npti, nptidx(1:npti), sst_1d(1:npti), sst_m(:,:) )
838         ! the following fields are modified in this routine
839         !!CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti), ato_i(:,:) )
840         !!CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d(1:npti,1:jpl), a_i(:,:,:) )
841         !!CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d  (1:npti,1:jpl), v_i  (:,:,:) )
842         CALL tab_3d_2d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s (:,:,:) )
843         CALL tab_3d_2d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
844         CALL tab_3d_2d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i(:,:,:) )
845         CALL tab_3d_2d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) )
846         CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) )
847         DO jl = 1, jpl
848            DO jk = 1, nlay_s
849               CALL tab_2d_1d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) )
850            END DO
851            DO jk = 1, nlay_i
852               CALL tab_2d_1d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
853            END DO
854         END DO
855         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_dyn_1d    (1:npti), sfx_dyn    (:,:) )
856         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_bri_1d    (1:npti), sfx_bri    (:,:) )
857         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_dyn_1d    (1:npti), wfx_dyn    (:,:) )
858         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_dyn_1d    (1:npti), hfx_dyn    (:,:) )
859         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) )
860         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d    (1:npti), wfx_pnd    (:,:) )
861         !
862         !                 !---------------------!
863      CASE( 2 )            !==  from 1D to 2D  ==!
864         !                 !---------------------!
865         CALL tab_1d_2d( npti, nptidx(1:npti), ato_i_1d(1:npti), ato_i(:,:) )
866         CALL tab_2d_3d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i (:,:,:) )
867         CALL tab_2d_3d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) )
868         CALL tab_2d_3d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s (:,:,:) )
869         CALL tab_2d_3d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
870         CALL tab_2d_3d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i(:,:,:) )
871         CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) )
872         CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) )
873         DO jl = 1, jpl
874            DO jk = 1, nlay_s
875               CALL tab_1d_2d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) )
876            END DO
877            DO jk = 1, nlay_i
878               CALL tab_1d_2d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
879            END DO
880         END DO
881         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_dyn_1d    (1:npti), sfx_dyn    (:,:) )
882         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_bri_1d    (1:npti), sfx_bri    (:,:) )
883         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_dyn_1d    (1:npti), wfx_dyn    (:,:) )
884         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_dyn_1d    (1:npti), hfx_dyn    (:,:) )
885         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) )
886         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d    (1:npti), wfx_pnd    (:,:) )
887         !
888      END SELECT
889      !
890   END SUBROUTINE ice_dyn_1d2d
891   
892
893   SUBROUTINE ice_dyn_rdgrft_init
894      !!-------------------------------------------------------------------
895      !!                  ***  ROUTINE ice_dyn_rdgrft_init ***
896      !!
897      !! ** Purpose :   Physical constants and parameters linked
898      !!                to the mechanical ice redistribution
899      !!
900      !! ** Method  :   Read the namdyn_rdgrft namelist
901      !!                and check the parameters values
902      !!                called at the first timestep (nit000)
903      !!
904      !! ** input   :   Namelist namdyn_rdgrft
905      !!-------------------------------------------------------------------
906      INTEGER :: ios                 ! Local integer output status for namelist read
907      !!
908      NAMELIST/namdyn_rdgrft/ ln_str_H79, rn_pstar, rn_crhg, &
909         &                    rn_csrdg  ,                    &
910         &                    ln_partf_lin, rn_gstar,        &
911         &                    ln_partf_exp, rn_astar,        & 
912         &                    ln_ridging, rn_hstar, rn_porordg, rn_fsnwrdg, rn_fpndrdg,  & 
913         &                    ln_rafting, rn_hraft, rn_craft  , rn_fsnwrft, rn_fpndrft
914      !!-------------------------------------------------------------------
915      !
916      REWIND( numnam_ice_ref )              ! Namelist namicetdme in reference namelist : Ice mechanical ice redistribution
917      READ  ( numnam_ice_ref, namdyn_rdgrft, IOSTAT = ios, ERR = 901)
918901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdyn_rdgrft in reference namelist', lwp )
919      REWIND( numnam_ice_cfg )              ! Namelist namdyn_rdgrft in configuration namelist : Ice mechanical ice redistribution
920      READ  ( numnam_ice_cfg, namdyn_rdgrft, IOSTAT = ios, ERR = 902 )
921902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdyn_rdgrft in configuration namelist', lwp )
922      IF(lwm) WRITE ( numoni, namdyn_rdgrft )
923      !
924      IF (lwp) THEN                          ! control print
925         WRITE(numout,*)
926         WRITE(numout,*) 'ice_dyn_rdgrft_init: ice parameters for ridging/rafting '
927         WRITE(numout,*) '~~~~~~~~~~~~~~~~~~'
928         WRITE(numout,*) '   Namelist namdyn_rdgrft:'
929         WRITE(numout,*) '      ice strength parameterization Hibler (1979)              ln_str_H79   = ', ln_str_H79 
930         WRITE(numout,*) '            1st bulk-rheology parameter                        rn_pstar     = ', rn_pstar
931         WRITE(numout,*) '            2nd bulk-rhelogy parameter                         rn_crhg      = ', rn_crhg
932         WRITE(numout,*) '      Fraction of shear energy contributing to ridging         rn_csrdg     = ', rn_csrdg 
933         WRITE(numout,*) '      linear ridging participation function                    ln_partf_lin = ', ln_partf_lin
934         WRITE(numout,*) '            Fraction of ice coverage contributing to ridging   rn_gstar     = ', rn_gstar
935         WRITE(numout,*) '      Exponential ridging participation function               ln_partf_exp = ', ln_partf_exp
936         WRITE(numout,*) '            Equivalent to G* for an exponential function       rn_astar     = ', rn_astar
937         WRITE(numout,*) '      Ridging of ice sheets or not                             ln_ridging   = ', ln_ridging
938         WRITE(numout,*) '            max ridged ice thickness                           rn_hstar     = ', rn_hstar
939         WRITE(numout,*) '            Initial porosity of ridges                         rn_porordg   = ', rn_porordg
940         WRITE(numout,*) '            Fraction of snow volume conserved during ridging   rn_fsnwrdg   = ', rn_fsnwrdg 
941         WRITE(numout,*) '            Fraction of pond volume conserved during ridging   rn_fpndrdg   = ', rn_fpndrdg 
942         WRITE(numout,*) '      Rafting of ice sheets or not                             ln_rafting   = ', ln_rafting
943         WRITE(numout,*) '            Parmeter thickness (threshold between ridge-raft)  rn_hraft     = ', rn_hraft
944         WRITE(numout,*) '            Rafting hyperbolic tangent coefficient             rn_craft     = ', rn_craft 
945         WRITE(numout,*) '            Fraction of snow volume conserved during rafting   rn_fsnwrft   = ', rn_fsnwrft 
946         WRITE(numout,*) '            Fraction of pond volume conserved during rafting   rn_fpndrft   = ', rn_fpndrft 
947      ENDIF
948      !
949      IF ( ( ln_partf_lin .AND. ln_partf_exp ) .OR. ( .NOT.ln_partf_lin .AND. .NOT.ln_partf_exp ) ) THEN
950         CALL ctl_stop( 'ice_dyn_rdgrft_init: choose one and only one participation function (ln_partf_lin or ln_partf_exp)' )
951      ENDIF
952      !
953      IF( .NOT. ln_icethd ) THEN
954         rn_porordg = 0._wp
955         rn_fsnwrdg = 1._wp ; rn_fsnwrft = 1._wp
956         rn_fpndrdg = 1._wp ; rn_fpndrft = 1._wp
957         IF( lwp ) THEN
958            WRITE(numout,*) '      ==> only ice dynamics is activated, thus some parameters must be changed'
959            WRITE(numout,*) '            rn_porordg   = ', rn_porordg
960            WRITE(numout,*) '            rn_fsnwrdg   = ', rn_fsnwrdg 
961            WRITE(numout,*) '            rn_fpndrdg   = ', rn_fpndrdg 
962            WRITE(numout,*) '            rn_fsnwrft   = ', rn_fsnwrft 
963            WRITE(numout,*) '            rn_fpndrft   = ', rn_fpndrft 
964         ENDIF
965      ENDIF
966      !                              ! allocate arrays
967      IF( ice_dyn_rdgrft_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'ice_dyn_rdgrft_init: unable to allocate arrays' )
968      !
969  END SUBROUTINE ice_dyn_rdgrft_init
970
971#else
972   !!----------------------------------------------------------------------
973   !!   Default option         Empty module           NO SI3 sea-ice model
974   !!----------------------------------------------------------------------
975#endif
976
977   !!======================================================================
978END MODULE icedyn_rdgrft
Note: See TracBrowser for help on using the repository browser.