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.
icedyn_rdgrft.F90 in NEMO/branches/2019/dev_r11842_SI3-10_EAP/src/ICE – NEMO

source: NEMO/branches/2019/dev_r11842_SI3-10_EAP/src/ICE/icedyn_rdgrft.F90 @ 12261

Last change on this file since 12261 was 12261, checked in by stefryn, 4 years ago

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