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/releases/r4.0/r4.0-HEAD/src/ICE – NEMO

source: NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icedyn_rdgrft.F90 @ 13494

Last change on this file since 13494 was 13494, checked in by clem, 4 years ago

4.0-HEAD: avoid one unecessary calculation, almost cosmetics
  • Property svn:keywords set to Id
File size: 52.5 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, zdelt  ! 1D divu_i & delta_i
140      !
141      INTEGER, PARAMETER ::   jp_itermax = 20   
142      !!-------------------------------------------------------------------
143      ! controls
144      IF( ln_timing    )   CALL timing_start('icedyn_rdgrft')                                                             ! timing
145      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icedyn_rdgrft', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
146      IF( ln_icediachk )   CALL ice_cons2D  (0, 'icedyn_rdgrft',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_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), zdelt   (1:npti)      , delta_i )
177         ! needed here and in the iteration loop
178         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)
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            IF( zdivu(ji) < 0._wp )   closing_net(ji) = MAX( closing_net(ji), -zdivu(ji) )   ! make sure the closing rate is large enough
189            !                                                                                ! to give asum = 1.0 after ridging
190            ! Opening rate (non-negative) that will give asum = 1.0 after ridging.
191            opning(ji) = closing_net(ji) + zdivu(ji)
192         END DO
193         !
194         !------------------------------------
195         ! 2) Identify grid cells with ridging
196         !------------------------------------
197         CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d, closing_net )
198
199         DO ji = 1, npti
200            IF( SUM( apartf(ji,1:jpl) ) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
201               ipti = ipti + 1
202               iptidx     (ipti)   = nptidx     (ji)
203               ! adjust to new indices
204               a_i_2d     (ipti,:) = a_i_2d     (ji,:)
205               v_i_2d     (ipti,:) = v_i_2d     (ji,:)
206               ato_i_1d   (ipti)   = ato_i_1d   (ji)
207               closing_net(ipti)   = closing_net(ji)
208               zdivu      (ipti)   = zdivu      (ji)
209               opning     (ipti)   = opning     (ji)
210            ENDIF
211         END DO
212
213      ENDIF
214
215      ! grid cells with ridging
216      nptidx(:) = iptidx(:)
217      npti      = ipti
218
219      !-----------------
220      ! 3) Start ridging
221      !-----------------
222      IF( npti > 0 ) THEN
223         
224         CALL ice_dyn_1d2d( 1 )            ! --- Move to 1D arrays --- !
225
226         iter            = 1
227         iterate_ridging = 1     
228         !                                                        !----------------------!
229         DO WHILE( iterate_ridging > 0 .AND. iter < jp_itermax )  !  ridging iterations  !
230            !                                                     !----------------------!
231            ! Calculate participation function (apartf)
232            !       and transfer      function
233            !       and closing_gross (+correction on opening)
234            CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d, closing_net )
235
236            ! Redistribute area, volume, and energy between categories
237            CALL rdgrft_shift
238
239            ! Do we keep on iterating?
240            !-------------------------
241            ! Check whether a_i + ato_i = 0
242            ! If not, because the closing and opening rates were reduced above, ridge again with new rates
243            iterate_ridging = 0
244            DO ji = 1, npti
245               zfac = 1._wp - ( ato_i_1d(ji) + SUM( a_i_2d(ji,:) ) )
246               IF( ABS( zfac ) < epsi10 ) THEN
247                  closing_net(ji) = 0._wp
248                  opning     (ji) = 0._wp
249                  ato_i_1d   (ji) = MAX( 0._wp, 1._wp - SUM( a_i_2d(ji,:) ) )
250               ELSE
251                  iterate_ridging  = 1
252                  zdivu      (ji) = zfac * r1_rdtice
253                  closing_net(ji) = MAX( 0._wp, -zdivu(ji) )
254                  opning     (ji) = MAX( 0._wp,  zdivu(ji) )
255               ENDIF
256            END DO
257            !
258            iter = iter + 1
259            IF( iter  >  jp_itermax )    CALL ctl_stop( 'STOP',  'icedyn_rdgrft: non-converging ridging scheme'  )
260            !
261         END DO
262
263         CALL ice_dyn_1d2d( 2 )            ! --- Move to 2D arrays --- !
264
265      ENDIF
266   
267      CALL ice_var_agg( 1 ) 
268
269      ! controls
270      IF( ln_ctl       )   CALL ice_prt3D   ('icedyn_rdgrft')                                                             ! prints
271      IF( ln_icectl    )   CALL ice_prt     (kt, iiceprt, jiceprt,-1, ' - ice dyn rdgrft - ')                             ! prints
272      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icedyn_rdgrft', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
273      IF( ln_icediachk )   CALL ice_cons2D  (1, 'icedyn_rdgrft',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) ! conservation
274      IF( ln_timing    )   CALL timing_stop ('icedyn_rdgrft')                                                             ! timing
275      !
276   END SUBROUTINE ice_dyn_rdgrft
277
278
279   SUBROUTINE rdgrft_prep( pa_i, pv_i, pato_i, pclosing_net )
280      !!-------------------------------------------------------------------
281      !!                ***  ROUTINE rdgrft_prep ***
282      !!
283      !! ** Purpose :   preparation for ridging calculations
284      !!
285      !! ** Method  :   Compute the thickness distribution of the ice and open water
286      !!                participating in ridging and of the resulting ridges.
287      !!-------------------------------------------------------------------
288      REAL(wp), DIMENSION(:)  , INTENT(in) ::   pato_i, pclosing_net 
289      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pa_i, pv_i 
290      !!
291      INTEGER  ::   ji, jl                     ! dummy loop indices
292      REAL(wp) ::   z1_gstar, z1_astar, zhmean, zfac   ! local scalar
293      REAL(wp), DIMENSION(jpij)        ::   zasum, z1_asum, zaksum   ! sum of a_i+ato_i and reverse
294      REAL(wp), DIMENSION(jpij,jpl)    ::   zhi                      ! ice thickness
295      REAL(wp), DIMENSION(jpij,-1:jpl) ::   zGsum                    ! zGsum(n) = sum of areas in categories 0 to n
296      !--------------------------------------------------------------------
297
298      z1_gstar = 1._wp / rn_gstar
299      z1_astar = 1._wp / rn_astar
300
301      !                       ! Ice thickness needed for rafting
302      WHERE( pa_i(1:npti,:) > epsi10 )   ;   zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:)
303      ELSEWHERE                          ;   zhi(1:npti,:) = 0._wp
304      END WHERE
305
306      ! 1) Participation function (apartf): a(h) = b(h).g(h)
307      !-----------------------------------------------------------------
308      ! Compute the participation function = total area lost due to ridging/closing
309      ! This is analogous to
310      !   a(h) = b(h)g(h) as defined in Thorndike et al. (1975).
311      !   assuming b(h) = (2/Gstar) * (1 - G(h)/Gstar).
312      !
313      ! apartf = integrating b(h)g(h) between the category boundaries
314      ! apartf is always >= 0 and SUM(apartf(0:jpl))=1
315      !-----------------------------------------------------------------
316      !
317      ! Compute total area of ice plus open water.
318      ! This is in general not equal to one because of divergence during transport
319      zasum(1:npti) = pato_i(1:npti) + SUM( pa_i(1:npti,:), dim=2 )
320      !
321      WHERE( zasum(1:npti) > epsi10 )   ;   z1_asum(1:npti) = 1._wp / zasum(1:npti)
322      ELSEWHERE                         ;   z1_asum(1:npti) = 0._wp
323      END WHERE
324      !
325      ! Compute cumulative thickness distribution function
326      ! Compute the cumulative thickness distribution function zGsum,
327      ! where zGsum(n) is the fractional area in categories 0 to n.
328      ! initial value (in h = 0) = open water area
329      zGsum(1:npti,-1) = 0._wp
330      zGsum(1:npti,0 ) = pato_i(1:npti) * z1_asum(1:npti)
331      DO jl = 1, jpl
332         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)
333      END DO
334      !
335      IF( ln_partf_lin ) THEN          !--- Linear formulation (Thorndike et al., 1975)
336         DO jl = 0, jpl   
337            DO ji = 1, npti
338               IF    ( zGsum(ji,jl)   < rn_gstar ) THEN
339                  apartf(ji,jl) = z1_gstar * ( zGsum(ji,jl) - zGsum(ji,jl-1) ) * &
340                     &                       ( 2._wp - ( zGsum(ji,jl-1) + zGsum(ji,jl) ) * z1_gstar )
341               ELSEIF( zGsum(ji,jl-1) < rn_gstar ) THEN
342                  apartf(ji,jl) = z1_gstar * ( rn_gstar     - zGsum(ji,jl-1) ) *  &
343                     &                       ( 2._wp - ( zGsum(ji,jl-1) + rn_gstar        ) * z1_gstar )
344               ELSE
345                  apartf(ji,jl) = 0._wp
346               ENDIF
347            END DO
348         END DO
349         !
350      ELSEIF( ln_partf_exp ) THEN      !--- Exponential, more stable formulation (Lipscomb et al, 2007)
351         !                       
352         zfac = 1._wp / ( 1._wp - EXP(-z1_astar) )
353         DO jl = -1, jpl
354            DO ji = 1, npti
355               zGsum(ji,jl) = EXP( -zGsum(ji,jl) * z1_astar ) * zfac
356            END DO
357         END DO
358         DO jl = 0, jpl
359            DO ji = 1, npti
360               apartf(ji,jl) = zGsum(ji,jl-1) - zGsum(ji,jl)
361            END DO
362         END DO
363         !
364      ENDIF
365
366      !                                !--- Ridging and rafting participation concentrations
367      IF( ln_rafting .AND. ln_ridging ) THEN             !- ridging & rafting
368         DO jl = 1, jpl
369            DO ji = 1, npti
370               aridge(ji,jl) = ( 1._wp + TANH ( rn_craft * ( zhi(ji,jl) - rn_hraft ) ) ) * 0.5_wp * apartf(ji,jl)
371               araft (ji,jl) = apartf(ji,jl) - aridge(ji,jl)
372            END DO
373         END DO
374      ELSEIF( ln_ridging .AND. .NOT. ln_rafting ) THEN   !- ridging alone
375         DO jl = 1, jpl
376            DO ji = 1, npti
377               aridge(ji,jl) = apartf(ji,jl)
378               araft (ji,jl) = 0._wp
379            END DO
380         END DO
381      ELSEIF( ln_rafting .AND. .NOT. ln_ridging ) THEN   !- rafting alone   
382         DO jl = 1, jpl
383            DO ji = 1, npti
384               aridge(ji,jl) = 0._wp
385               araft (ji,jl) = apartf(ji,jl)
386            END DO
387         END DO
388      ELSE                                               !- no ridging & no rafting
389         DO jl = 1, jpl
390            DO ji = 1, npti
391               aridge(ji,jl) = 0._wp
392               araft (ji,jl) = 0._wp         
393            END DO
394         END DO
395      ENDIF
396
397      ! 2) Transfer function
398      !-----------------------------------------------------------------
399      ! Compute max and min ridged ice thickness for each ridging category.
400      ! Assume ridged ice is uniformly distributed between hrmin and hrmax.
401      !
402      ! This parameterization is a modified version of Hibler (1980).
403      ! The mean ridging thickness, zhmean, is proportional to hi^(0.5)
404      !  and for very thick ridging ice must be >= hrdg_hi_min*hi
405      !
406      ! The minimum ridging thickness, hrmin, is equal to 2*hi
407      !  (i.e., rafting) and for very thick ridging ice is
408      !  constrained by hrmin <= (zhmean + hi)/2.
409      !
410      ! The maximum ridging thickness, hrmax, is determined by zhmean and hrmin.
411      !
412      ! These modifications have the effect of reducing the ice strength
413      ! (relative to the Hibler formulation) when very thick ice is ridging.
414      !
415      ! zaksum = net area removed/ total area removed
416      ! where total area removed = area of ice that ridges
417      !         net area removed = total area removed - area of new ridges
418      !-----------------------------------------------------------------
419      zfac = 1._wp / hi_hrft
420      zaksum(1:npti) = apartf(1:npti,0)
421      !
422      DO jl = 1, jpl
423         DO ji = 1, npti
424            IF ( apartf(ji,jl) > 0._wp ) THEN
425               zhmean         = MAX( SQRT( rn_hstar * zhi(ji,jl) ), zhi(ji,jl) * hrdg_hi_min )
426               hrmin  (ji,jl) = MIN( 2._wp * zhi(ji,jl), 0.5_wp * ( zhmean + zhi(ji,jl) ) )
427               hrmax  (ji,jl) = 2._wp * zhmean - hrmin(ji,jl)
428               hraft  (ji,jl) = zhi(ji,jl) * zfac
429               hi_hrdg(ji,jl) = zhi(ji,jl) / MAX( zhmean, epsi20 )
430               !
431               ! Normalization factor : zaksum, ensures mass conservation
432               zaksum(ji) = zaksum(ji) + aridge(ji,jl) * ( 1._wp - hi_hrdg(ji,jl) )    &
433                  &                    + araft (ji,jl) * ( 1._wp - hi_hrft )
434            ELSE
435               hrmin  (ji,jl) = 0._wp 
436               hrmax  (ji,jl) = 0._wp 
437               hraft  (ji,jl) = 0._wp 
438               hi_hrdg(ji,jl) = 1._wp
439            ENDIF
440         END DO
441      END DO
442      !
443      ! 3) closing_gross
444      !-----------------
445      ! Based on the ITD of ridging and ridged ice, convert the net closing rate to a gross closing rate. 
446      ! NOTE: 0 < aksum <= 1
447      WHERE( zaksum(1:npti) > epsi10 )   ;   closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti)
448      ELSEWHERE                          ;   closing_gross(1:npti) = 0._wp
449      END WHERE
450     
451      ! correction to closing rate if excessive ice removal
452      !----------------------------------------------------
453      ! Reduce the closing rate if more than 100% of any ice category would be removed
454      ! Reduce the opening rate in proportion
455      DO jl = 1, jpl
456         DO ji = 1, npti
457            zfac = apartf(ji,jl) * closing_gross(ji) * rdt_ice
458            IF( zfac > pa_i(ji,jl) .AND. apartf(ji,jl) /= 0._wp ) THEN
459               closing_gross(ji) = pa_i(ji,jl) / apartf(ji,jl) * r1_rdtice
460            ENDIF
461         END DO
462      END DO     
463
464      ! 4) correction to opening if excessive open water removal
465      !---------------------------------------------------------
466      ! Reduce the closing rate if more than 100% of the open water would be removed
467      ! Reduce the opening rate in proportion
468      DO ji = 1, npti 
469         zfac = pato_i(ji) + ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * rdt_ice
470         IF( zfac < 0._wp ) THEN           ! would lead to negative ato_i
471            opning(ji) = apartf(ji,0) * closing_gross(ji) - pato_i(ji) * r1_rdtice 
472         ELSEIF( zfac > zasum(ji) ) THEN   ! would lead to ato_i > asum
473            opning(ji) = apartf(ji,0) * closing_gross(ji) + ( zasum(ji) - pato_i(ji) ) * r1_rdtice 
474         ENDIF
475      END DO
476      !
477   END SUBROUTINE rdgrft_prep
478
479
480   SUBROUTINE rdgrft_shift
481      !!-------------------------------------------------------------------
482      !!                ***  ROUTINE rdgrft_shift ***
483      !!
484      !! ** Purpose :   shift ridging ice among thickness categories of ice thickness
485      !!
486      !! ** Method  :   Remove area, volume, and energy from each ridging category
487      !!                and add to thicker ice categories.
488      !!-------------------------------------------------------------------
489      !
490      INTEGER  ::   ji, jj, jl, jl1, jl2, jk   ! dummy loop indices
491      REAL(wp) ::   hL, hR, farea              ! left and right limits of integration and new area going to jl2
492      REAL(wp) ::   vsw                        ! vol of water trapped into ridges
493      REAL(wp) ::   afrdg, afrft               ! fraction of category area ridged/rafted
494      REAL(wp)                  ::   airdg1, oirdg1, aprdg1, virdg1, sirdg1
495      REAL(wp)                  ::   airft1, oirft1, aprft1
496      REAL(wp), DIMENSION(jpij) ::   airdg2, oirdg2, aprdg2, virdg2, sirdg2, vsrdg, vprdg, vlrdg  ! area etc of new ridges
497      REAL(wp), DIMENSION(jpij) ::   airft2, oirft2, aprft2, virft , sirft , vsrft, vprft, vlrft  ! area etc of rafted ice
498      !
499      REAL(wp), DIMENSION(jpij) ::   ersw             ! enth of water trapped into ridges
500      REAL(wp), DIMENSION(jpij) ::   zswitch, fvol    ! new ridge volume going to jl2
501      REAL(wp), DIMENSION(jpij) ::   z1_ai            ! 1 / a
502      REAL(wp), DIMENSION(jpij) ::   zvti             ! sum(v_i)
503      !
504      REAL(wp), DIMENSION(jpij,nlay_s) ::   esrft     ! snow energy of rafting ice
505      REAL(wp), DIMENSION(jpij,nlay_i) ::   eirft     ! ice  energy of rafting ice
506      REAL(wp), DIMENSION(jpij,nlay_s) ::   esrdg     ! enth*volume of new ridges     
507      REAL(wp), DIMENSION(jpij,nlay_i) ::   eirdg     ! enth*volume of new ridges
508      !
509      INTEGER , DIMENSION(jpij) ::   itest_rdg, itest_rft   ! test for conservation
510      !!-------------------------------------------------------------------
511      !
512      zvti(1:npti) = SUM( v_i_2d(1:npti,:), dim=2 )   ! total ice volume
513      !
514      ! 1) Change in open water area due to closing and opening
515      !--------------------------------------------------------
516      DO ji = 1, npti
517         ato_i_1d(ji) = MAX( 0._wp, ato_i_1d(ji) + ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * rdt_ice )
518      END DO
519     
520      ! 2) compute categories in which ice is removed (jl1)
521      !----------------------------------------------------
522      DO jl1 = 1, jpl
523
524         IF( nn_icesal /= 2 )  THEN     
525            CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d(1:npti), s_i(:,:,jl1) )
526         ENDIF
527
528         DO ji = 1, npti
529
530            IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN   ! only if ice is ridging
531
532               IF( a_i_2d(ji,jl1) > epsi10 ) THEN   ;   z1_ai(ji) = 1._wp / a_i_2d(ji,jl1)
533               ELSE                                 ;   z1_ai(ji) = 0._wp
534               ENDIF
535               
536               ! area of ridging / rafting ice (airdg1) and of new ridge (airdg2)
537               airdg1 = aridge(ji,jl1) * closing_gross(ji) * rdt_ice
538               airft1 = araft (ji,jl1) * closing_gross(ji) * rdt_ice
539
540               airdg2(ji) = airdg1 * hi_hrdg(ji,jl1)
541               airft2(ji) = airft1 * hi_hrft
542
543               ! ridging /rafting fractions
544               afrdg = airdg1 * z1_ai(ji)
545               afrft = airft1 * z1_ai(ji)
546
547               ! volume and enthalpy (J/m2, >0) of seawater trapped into ridges
548               IF    ( zvti(ji) <= 10. ) THEN ; vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg                                           ! v <= 10m then porosity = rn_porordg
549               ELSEIF( zvti(ji) >= 20. ) THEN ; vsw = 0._wp                                                                         ! v >= 20m then porosity = 0
550               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
551               ENDIF
552               ersw(ji) = -rhoi * vsw * rcp * sst_1d(ji)   ! clem: if sst>0, then ersw <0 (is that possible?)
553
554               ! volume etc of ridging / rafting ice and new ridges (vi, vs, sm, oi, es, ei)
555               virdg1     = v_i_2d (ji,jl1)   * afrdg
556               virdg2(ji) = v_i_2d (ji,jl1)   * afrdg + vsw
557               vsrdg(ji)  = v_s_2d (ji,jl1)   * afrdg
558               sirdg1     = sv_i_2d(ji,jl1)   * afrdg
559               sirdg2(ji) = sv_i_2d(ji,jl1)   * afrdg + vsw * sss_1d(ji)
560               oirdg1     = oa_i_2d(ji,jl1)   * afrdg
561               oirdg2(ji) = oa_i_2d(ji,jl1)   * afrdg * hi_hrdg(ji,jl1) 
562
563               virft(ji)  = v_i_2d (ji,jl1)   * afrft
564               vsrft(ji)  = v_s_2d (ji,jl1)   * afrft
565               sirft(ji)  = sv_i_2d(ji,jl1)   * afrft 
566               oirft1     = oa_i_2d(ji,jl1)   * afrft 
567               oirft2(ji) = oa_i_2d(ji,jl1)   * afrft * hi_hrft 
568
569               IF ( ln_pnd_LEV ) THEN
570                  aprdg1     = a_ip_2d(ji,jl1) * afrdg
571                  aprdg2(ji) = a_ip_2d(ji,jl1) * afrdg * hi_hrdg(ji,jl1)
572                  vprdg (ji) = v_ip_2d(ji,jl1) * afrdg
573                  aprft1     = a_ip_2d(ji,jl1) * afrft
574                  aprft2(ji) = a_ip_2d(ji,jl1) * afrft * hi_hrft
575                  vprft (ji) = v_ip_2d(ji,jl1) * afrft
576                  IF ( ln_pnd_lids ) THEN
577                     vlrdg (ji) = v_il_2d(ji,jl1) * afrdg
578                     vlrft (ji) = v_il_2d(ji,jl1) * afrft
579                  ENDIF
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_LEV ) 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                  IF ( ln_pnd_lids ) THEN
610                     v_il_2d(ji,jl1) = v_il_2d(ji,jl1) - vlrdg(ji) - vlrft(ji)
611                  ENDIF
612               ENDIF
613            ENDIF
614
615         END DO ! ji
616
617         ! special loop for e_s because of layers jk
618         DO jk = 1, nlay_s
619            DO ji = 1, npti
620               IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
621                  ! Compute ridging /rafting fractions
622                  afrdg = aridge(ji,jl1) * closing_gross(ji) * rdt_ice * z1_ai(ji)
623                  afrft = araft (ji,jl1) * closing_gross(ji) * rdt_ice * z1_ai(ji)
624                  ! Compute ridging /rafting ice and new ridges for es
625                  esrdg(ji,jk) = ze_s_2d (ji,jk,jl1) * afrdg
626                  esrft(ji,jk) = ze_s_2d (ji,jk,jl1) * afrft
627                  ! Put the snow lost by ridging into the ocean
628                  hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ( - esrdg(ji,jk) * ( 1._wp - rn_fsnwrdg )   &                 ! heat sink for ocean (<0, W.m-2)
629                     &                                - esrft(ji,jk) * ( 1._wp - rn_fsnwrft ) ) * r1_rdtice
630                  !
631                  ! Remove energy of new ridge to each category jl1
632                  !-------------------------------------------------
633                  ze_s_2d(ji,jk,jl1) = ze_s_2d(ji,jk,jl1) * ( 1._wp - afrdg - afrft ) 
634               ENDIF
635            END DO
636         END DO
637                 
638         ! special loop for e_i because of layers jk
639         DO jk = 1, nlay_i
640            DO ji = 1, npti
641               IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
642                  ! Compute ridging /rafting fractions
643                  afrdg = aridge(ji,jl1) * closing_gross(ji) * rdt_ice * z1_ai(ji)
644                  afrft = araft (ji,jl1) * closing_gross(ji) * rdt_ice * z1_ai(ji)
645                  ! Compute ridging ice and new ridges for ei
646                  eirdg(ji,jk) = ze_i_2d (ji,jk,jl1) * afrdg + ersw(ji) * r1_nlay_i
647                  eirft(ji,jk) = ze_i_2d (ji,jk,jl1) * afrft
648                  !
649                  ! Remove energy of new ridge to each category jl1
650                  !-------------------------------------------------
651                  ze_i_2d(ji,jk,jl1) = ze_i_2d(ji,jk,jl1) * ( 1._wp - afrdg - afrft ) 
652               ENDIF
653            END DO
654         END DO
655         
656         ! 3) compute categories in which ice is added (jl2)
657         !--------------------------------------------------
658         itest_rdg(1:npti) = 0
659         itest_rft(1:npti) = 0
660         DO jl2  = 1, jpl 
661            !
662            DO ji = 1, npti
663
664               IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
665
666                  ! Compute the fraction of ridged ice area and volume going to thickness category jl2
667                  IF( hrmin(ji,jl1) <= hi_max(jl2) .AND. hrmax(ji,jl1) > hi_max(jl2-1) ) THEN
668                     hL = MAX( hrmin(ji,jl1), hi_max(jl2-1) )
669                     hR = MIN( hrmax(ji,jl1), hi_max(jl2)   )
670                     farea    = ( hR      - hL      ) / ( hrmax(ji,jl1)                 - hrmin(ji,jl1)                 )
671                     fvol(ji) = ( hR * hR - hL * hL ) / ( hrmax(ji,jl1) * hrmax(ji,jl1) - hrmin(ji,jl1) * hrmin(ji,jl1) )
672                     !
673                     itest_rdg(ji) = 1   ! test for conservation
674                  ELSE
675                     farea    = 0._wp 
676                     fvol(ji) = 0._wp                 
677                  ENDIF
678
679                  ! Compute the fraction of rafted ice area and volume going to thickness category jl2
680                  IF( hraft(ji,jl1) <= hi_max(jl2) .AND. hraft(ji,jl1) >  hi_max(jl2-1) ) THEN
681                     zswitch(ji) = 1._wp
682                     !
683                     itest_rft(ji) = 1   ! test for conservation
684                  ELSE
685                     zswitch(ji) = 0._wp
686                  ENDIF
687                  !
688                  ! Patch to ensure perfect conservation if ice thickness goes mad
689                  ! Sometimes thickness is larger than hi_max(jpl) because of advection scheme (for very small areas)
690                  ! Then ice volume is removed from one category but the ridging/rafting scheme
691                  ! does not know where to move it, leading to a conservation issue. 
692                  IF( itest_rdg(ji) == 0 .AND. jl2 == jpl ) THEN   ;   farea = 1._wp   ;   fvol(ji) = 1._wp   ;   ENDIF
693                  IF( itest_rft(ji) == 0 .AND. jl2 == jpl )      zswitch(ji) = 1._wp
694                  !
695                  ! Add area, volume of new ridge to category jl2
696                  !----------------------------------------------
697                  a_i_2d (ji,jl2) = a_i_2d (ji,jl2) + ( airdg2(ji) * farea    + airft2(ji) * zswitch(ji) )
698                  oa_i_2d(ji,jl2) = oa_i_2d(ji,jl2) + ( oirdg2(ji) * farea    + oirft2(ji) * zswitch(ji) )
699                  v_i_2d (ji,jl2) = v_i_2d (ji,jl2) + ( virdg2(ji) * fvol(ji) + virft (ji) * zswitch(ji) )
700                  sv_i_2d(ji,jl2) = sv_i_2d(ji,jl2) + ( sirdg2(ji) * fvol(ji) + sirft (ji) * zswitch(ji) )
701                  v_s_2d (ji,jl2) = v_s_2d (ji,jl2) + ( vsrdg (ji) * rn_fsnwrdg * fvol(ji)  +  &
702                     &                                  vsrft (ji) * rn_fsnwrft * zswitch(ji) )
703                  IF ( ln_pnd_LEV ) THEN
704                     v_ip_2d (ji,jl2) = v_ip_2d(ji,jl2) + (   vprdg (ji) * rn_fpndrdg * fvol   (ji)   &
705                        &                                   + vprft (ji) * rn_fpndrft * zswitch(ji)   )
706                     a_ip_2d (ji,jl2) = a_ip_2d(ji,jl2) + (   aprdg2(ji) * rn_fpndrdg * farea         & 
707                        &                                   + aprft2(ji) * rn_fpndrft * zswitch(ji)   )
708                     IF ( ln_pnd_lids ) THEN
709                        v_il_2d (ji,jl2) = v_il_2d(ji,jl2) + (   vlrdg(ji) * rn_fpndrdg * fvol   (ji) &
710                           &                                   + vlrft(ji) * rn_fpndrft * zswitch(ji) )
711                     ENDIF
712                  ENDIF
713                 
714               ENDIF
715
716            END DO
717            ! Add snow energy of new ridge to category jl2
718            !---------------------------------------------
719            DO jk = 1, nlay_s
720               DO ji = 1, npti
721                  IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp )   &
722                     &   ze_s_2d(ji,jk,jl2) = ze_s_2d(ji,jk,jl2) + ( esrdg(ji,jk) * rn_fsnwrdg * fvol(ji)  +  &
723                     &                                               esrft(ji,jk) * rn_fsnwrft * zswitch(ji) )
724               END DO
725            END DO
726            ! Add ice energy of new ridge to category jl2
727            !--------------------------------------------
728            DO jk = 1, nlay_i
729               DO ji = 1, npti
730                  IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp )   &
731                     &   ze_i_2d(ji,jk,jl2) = ze_i_2d(ji,jk,jl2) + eirdg(ji,jk) * fvol(ji) + eirft(ji,jk) * zswitch(ji)                 
732               END DO
733            END DO
734            !
735         END DO ! jl2
736         !
737      END DO ! jl1
738      !
739      ! roundoff errors
740      !----------------
741      ! In case ridging/rafting lead to very small negative values (sometimes it happens)
742      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, v_il_2d, ze_s_2d, ze_i_2d )
743      !
744   END SUBROUTINE rdgrft_shift
745
746
747   SUBROUTINE ice_strength
748      !!----------------------------------------------------------------------
749      !!                ***  ROUTINE ice_strength ***
750      !!
751      !! ** Purpose :   computes ice strength used in dynamics routines of ice thickness
752      !!
753      !! ** Method  :   Compute the strength of the ice pack, defined as the energy (J m-2)
754      !!              dissipated per unit area removed from the ice pack under compression,
755      !!              and assumed proportional to the change in potential energy caused
756      !!              by ridging. Note that only Hibler's formulation is stable and that
757      !!              ice strength has to be smoothed
758      !!----------------------------------------------------------------------
759      INTEGER             ::   ji, jj, jl  ! dummy loop indices
760      INTEGER             ::   ismooth     ! smoothing the resistance to deformation
761      INTEGER             ::   itframe     ! number of time steps for the P smoothing
762      REAL(wp)            ::   zp, z1_3    ! local scalars
763      REAL(wp), DIMENSION(jpi,jpj) ::   zworka           ! temporary array used here
764      REAL(wp), DIMENSION(jpi,jpj) ::   zstrp1, zstrp2   ! strength at previous time steps
765      !!----------------------------------------------------------------------
766      !                              !--------------------------------------------------!
767      IF( ln_str_H79 ) THEN          ! Ice strength => Hibler (1979) method             !
768      !                              !--------------------------------------------------!
769         strength(:,:) = rn_pstar * SUM( v_i(:,:,:), dim=3 ) * EXP( -rn_crhg * ( 1._wp - SUM( a_i(:,:,:), dim=3 ) ) )
770         ismooth = 1
771         !                           !--------------------------------------------------!
772      ELSE                           ! Zero strength                                    !
773         !                           !--------------------------------------------------!
774         strength(:,:) = 0._wp
775         ismooth = 0
776      ENDIF
777      !                              !--------------------------------------------------!
778      SELECT CASE( ismooth )         ! Smoothing ice strength                           !
779      !                              !--------------------------------------------------!
780      CASE( 1 )               !--- Spatial smoothing
781         DO jj = 2, jpjm1
782            DO ji = 2, jpim1
783               IF ( SUM( a_i(ji,jj,:) ) > 0._wp ) THEN
784                  zworka(ji,jj) = ( 4.0 * strength(ji,jj)              &
785                     &                  + strength(ji-1,jj) * tmask(ji-1,jj,1) + strength(ji+1,jj) * tmask(ji+1,jj,1) & 
786                     &                  + strength(ji,jj-1) * tmask(ji,jj-1,1) + strength(ji,jj+1) * tmask(ji,jj+1,1) &
787                     &            ) / ( 4.0 + tmask(ji-1,jj,1) + tmask(ji+1,jj,1) + tmask(ji,jj-1,1) + tmask(ji,jj+1,1) )
788               ELSE
789                  zworka(ji,jj) = 0._wp
790               ENDIF
791            END DO
792         END DO
793         
794         DO jj = 2, jpjm1
795            DO ji = 2, jpim1
796               strength(ji,jj) = zworka(ji,jj)
797            END DO
798         END DO
799         CALL lbc_lnk( 'icedyn_rdgrft', strength, 'T', 1. )
800         !
801      CASE( 2 )               !--- Temporal smoothing
802         IF ( kt_ice == nit000 ) THEN
803            zstrp1(:,:) = 0._wp
804            zstrp2(:,:) = 0._wp
805         ENDIF
806         !
807         DO jj = 2, jpjm1
808            DO ji = 2, jpim1
809               IF ( SUM( a_i(ji,jj,:) ) > 0._wp ) THEN
810                  itframe = 1 ! number of time steps for the running mean
811                  IF ( zstrp1(ji,jj) > 0._wp ) itframe = itframe + 1
812                  IF ( zstrp2(ji,jj) > 0._wp ) itframe = itframe + 1
813                  zp = ( strength(ji,jj) + zstrp1(ji,jj) + zstrp2(ji,jj) ) / itframe
814                  zstrp2  (ji,jj) = zstrp1  (ji,jj)
815                  zstrp1  (ji,jj) = strength(ji,jj)
816                  strength(ji,jj) = zp
817               ENDIF
818            END DO
819         END DO
820         CALL lbc_lnk( 'icedyn_rdgrft', strength, 'T', 1. )
821         !
822      END SELECT
823      !
824   END SUBROUTINE ice_strength
825
826   
827   SUBROUTINE ice_dyn_1d2d( kn )
828      !!-----------------------------------------------------------------------
829      !!                   ***  ROUTINE ice_dyn_1d2d ***
830      !!                 
831      !! ** Purpose :   move arrays from 1d to 2d and the reverse
832      !!-----------------------------------------------------------------------
833      INTEGER, INTENT(in) ::   kn   ! 1= from 2D to 1D   ;   2= from 1D to 2D
834      !
835      INTEGER ::   jl, jk   ! dummy loop indices
836      !!-----------------------------------------------------------------------
837      !
838      SELECT CASE( kn )
839      !                    !---------------------!
840      CASE( 1 )            !==  from 2D to 1D  ==!
841         !                 !---------------------!
842         ! fields used but not modified
843         CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d(1:npti), sss_m(:,:) )
844         CALL tab_2d_1d( npti, nptidx(1:npti), sst_1d(1:npti), sst_m(:,:) )
845         ! the following fields are modified in this routine
846         !!CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti), ato_i(:,:) )
847         !!CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d(1:npti,1:jpl), a_i(:,:,:) )
848         !!CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d  (1:npti,1:jpl), v_i  (:,:,:) )
849         CALL tab_3d_2d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s (:,:,:) )
850         CALL tab_3d_2d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
851         CALL tab_3d_2d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i(:,:,:) )
852         CALL tab_3d_2d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) )
853         CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) )
854         CALL tab_3d_2d( npti, nptidx(1:npti), v_il_2d(1:npti,1:jpl), v_il(:,:,:) )
855         DO jl = 1, jpl
856            DO jk = 1, nlay_s
857               CALL tab_2d_1d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) )
858            END DO
859            DO jk = 1, nlay_i
860               CALL tab_2d_1d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
861            END DO
862         END DO
863         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_dyn_1d    (1:npti), sfx_dyn    (:,:) )
864         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_bri_1d    (1:npti), sfx_bri    (:,:) )
865         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_dyn_1d    (1:npti), wfx_dyn    (:,:) )
866         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_dyn_1d    (1:npti), hfx_dyn    (:,:) )
867         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) )
868         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d    (1:npti), wfx_pnd    (:,:) )
869         !
870         !                 !---------------------!
871      CASE( 2 )            !==  from 1D to 2D  ==!
872         !                 !---------------------!
873         CALL tab_1d_2d( npti, nptidx(1:npti), ato_i_1d(1:npti), ato_i(:,:) )
874         CALL tab_2d_3d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i (:,:,:) )
875         CALL tab_2d_3d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) )
876         CALL tab_2d_3d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s (:,:,:) )
877         CALL tab_2d_3d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
878         CALL tab_2d_3d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i(:,:,:) )
879         CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) )
880         CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) )
881         CALL tab_2d_3d( npti, nptidx(1:npti), v_il_2d(1:npti,1:jpl), v_il(:,:,:) )
882         DO jl = 1, jpl
883            DO jk = 1, nlay_s
884               CALL tab_1d_2d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) )
885            END DO
886            DO jk = 1, nlay_i
887               CALL tab_1d_2d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
888            END DO
889         END DO
890         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_dyn_1d    (1:npti), sfx_dyn    (:,:) )
891         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_bri_1d    (1:npti), sfx_bri    (:,:) )
892         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_dyn_1d    (1:npti), wfx_dyn    (:,:) )
893         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_dyn_1d    (1:npti), hfx_dyn    (:,:) )
894         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) )
895         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d    (1:npti), wfx_pnd    (:,:) )
896         !
897      END SELECT
898      !
899   END SUBROUTINE ice_dyn_1d2d
900   
901
902   SUBROUTINE ice_dyn_rdgrft_init
903      !!-------------------------------------------------------------------
904      !!                  ***  ROUTINE ice_dyn_rdgrft_init ***
905      !!
906      !! ** Purpose :   Physical constants and parameters linked
907      !!                to the mechanical ice redistribution
908      !!
909      !! ** Method  :   Read the namdyn_rdgrft namelist
910      !!                and check the parameters values
911      !!                called at the first timestep (nit000)
912      !!
913      !! ** input   :   Namelist namdyn_rdgrft
914      !!-------------------------------------------------------------------
915      INTEGER :: ios                 ! Local integer output status for namelist read
916      !!
917      NAMELIST/namdyn_rdgrft/ ln_str_H79, rn_pstar, rn_crhg, &
918         &                    rn_csrdg  ,                    &
919         &                    ln_partf_lin, rn_gstar,        &
920         &                    ln_partf_exp, rn_astar,        & 
921         &                    ln_ridging, rn_hstar, rn_porordg, rn_fsnwrdg, rn_fpndrdg,  & 
922         &                    ln_rafting, rn_hraft, rn_craft  , rn_fsnwrft, rn_fpndrft
923      !!-------------------------------------------------------------------
924      !
925      REWIND( numnam_ice_ref )              ! Namelist namicetdme in reference namelist : Ice mechanical ice redistribution
926      READ  ( numnam_ice_ref, namdyn_rdgrft, IOSTAT = ios, ERR = 901)
927901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdyn_rdgrft in reference namelist' )
928      REWIND( numnam_ice_cfg )              ! Namelist namdyn_rdgrft in configuration namelist : Ice mechanical ice redistribution
929      READ  ( numnam_ice_cfg, namdyn_rdgrft, IOSTAT = ios, ERR = 902 )
930902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdyn_rdgrft in configuration namelist' )
931      IF(lwm) WRITE ( numoni, namdyn_rdgrft )
932      !
933      IF (lwp) THEN                          ! control print
934         WRITE(numout,*)
935         WRITE(numout,*) 'ice_dyn_rdgrft_init: ice parameters for ridging/rafting '
936         WRITE(numout,*) '~~~~~~~~~~~~~~~~~~'
937         WRITE(numout,*) '   Namelist namdyn_rdgrft:'
938         WRITE(numout,*) '      ice strength parameterization Hibler (1979)              ln_str_H79   = ', ln_str_H79 
939         WRITE(numout,*) '            1st bulk-rheology parameter                        rn_pstar     = ', rn_pstar
940         WRITE(numout,*) '            2nd bulk-rhelogy parameter                         rn_crhg      = ', rn_crhg
941         WRITE(numout,*) '      Fraction of shear energy contributing to ridging         rn_csrdg     = ', rn_csrdg 
942         WRITE(numout,*) '      linear ridging participation function                    ln_partf_lin = ', ln_partf_lin
943         WRITE(numout,*) '            Fraction of ice coverage contributing to ridging   rn_gstar     = ', rn_gstar
944         WRITE(numout,*) '      Exponential ridging participation function               ln_partf_exp = ', ln_partf_exp
945         WRITE(numout,*) '            Equivalent to G* for an exponential function       rn_astar     = ', rn_astar
946         WRITE(numout,*) '      Ridging of ice sheets or not                             ln_ridging   = ', ln_ridging
947         WRITE(numout,*) '            max ridged ice thickness                           rn_hstar     = ', rn_hstar
948         WRITE(numout,*) '            Initial porosity of ridges                         rn_porordg   = ', rn_porordg
949         WRITE(numout,*) '            Fraction of snow volume conserved during ridging   rn_fsnwrdg   = ', rn_fsnwrdg 
950         WRITE(numout,*) '            Fraction of pond volume conserved during ridging   rn_fpndrdg   = ', rn_fpndrdg 
951         WRITE(numout,*) '      Rafting of ice sheets or not                             ln_rafting   = ', ln_rafting
952         WRITE(numout,*) '            Parmeter thickness (threshold between ridge-raft)  rn_hraft     = ', rn_hraft
953         WRITE(numout,*) '            Rafting hyperbolic tangent coefficient             rn_craft     = ', rn_craft 
954         WRITE(numout,*) '            Fraction of snow volume conserved during rafting   rn_fsnwrft   = ', rn_fsnwrft 
955         WRITE(numout,*) '            Fraction of pond volume conserved during rafting   rn_fpndrft   = ', rn_fpndrft 
956      ENDIF
957      !
958      IF ( ( ln_partf_lin .AND. ln_partf_exp ) .OR. ( .NOT.ln_partf_lin .AND. .NOT.ln_partf_exp ) ) THEN
959         CALL ctl_stop( 'ice_dyn_rdgrft_init: choose one and only one participation function (ln_partf_lin or ln_partf_exp)' )
960      ENDIF
961      !
962      IF( .NOT. ln_icethd ) THEN
963         rn_porordg = 0._wp
964         rn_fsnwrdg = 1._wp ; rn_fsnwrft = 1._wp
965         rn_fpndrdg = 1._wp ; rn_fpndrft = 1._wp
966         IF( lwp ) THEN
967            WRITE(numout,*) '      ==> only ice dynamics is activated, thus some parameters must be changed'
968            WRITE(numout,*) '            rn_porordg   = ', rn_porordg
969            WRITE(numout,*) '            rn_fsnwrdg   = ', rn_fsnwrdg 
970            WRITE(numout,*) '            rn_fpndrdg   = ', rn_fpndrdg 
971            WRITE(numout,*) '            rn_fsnwrft   = ', rn_fsnwrft 
972            WRITE(numout,*) '            rn_fpndrft   = ', rn_fpndrft 
973         ENDIF
974      ENDIF
975      !                              ! allocate arrays
976      IF( ice_dyn_rdgrft_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'ice_dyn_rdgrft_init: unable to allocate arrays' )
977      !
978  END SUBROUTINE ice_dyn_rdgrft_init
979
980#else
981   !!----------------------------------------------------------------------
982   !!   Default option         Empty module           NO SI3 sea-ice model
983   !!----------------------------------------------------------------------
984#endif
985
986   !!======================================================================
987END MODULE icedyn_rdgrft
Note: See TracBrowser for help on using the repository browser.