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/trunk/src/ICE – NEMO

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

Last change on this file since 13495 was 13495, checked in by clem, 20 months ago

trunk: 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   !! * Substitutions
78#  include "do_loop_substitute.h90"
79   !!----------------------------------------------------------------------
80   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
81   !! $Id$
82   !! Software governed by the CeCILL licence     (./LICENSE)
83   !!----------------------------------------------------------------------
84CONTAINS
85
86   INTEGER FUNCTION ice_dyn_rdgrft_alloc()
87      !!-------------------------------------------------------------------
88      !!                ***  ROUTINE ice_dyn_rdgrft_alloc ***
89      !!-------------------------------------------------------------------
90      ALLOCATE( closing_net(jpij)  , opning(jpij)      , closing_gross(jpij) ,               &
91         &      apartf(jpij,0:jpl) , hrmin  (jpij,jpl) , hraft(jpij,jpl) , aridge(jpij,jpl), &
92         &      hrmax (jpij,jpl)   , hi_hrdg(jpij,jpl) , araft(jpij,jpl) ,                   &
93         &      ze_i_2d(jpij,nlay_i,jpl), ze_s_2d(jpij,nlay_s,jpl), STAT=ice_dyn_rdgrft_alloc )
94
95      CALL mpp_sum ( 'icedyn_rdgrft', ice_dyn_rdgrft_alloc )
96      IF( ice_dyn_rdgrft_alloc /= 0 )   CALL ctl_stop( 'STOP',  'ice_dyn_rdgrft_alloc: failed to allocate arrays'  )
97      !
98   END FUNCTION ice_dyn_rdgrft_alloc
99
100
101   SUBROUTINE ice_dyn_rdgrft( kt )
102      !!-------------------------------------------------------------------
103      !!                ***  ROUTINE ice_dyn_rdgrft ***
104      !!
105      !! ** Purpose :   computes the mechanical redistribution of ice thickness
106      !!
107      !! ** Method  :   Steps :
108      !!       0) Identify grid cells with ice
109      !!       1) Calculate closing rate, divergence and opening
110      !!       2) Identify grid cells with ridging
111      !!       3) Start ridging iterations
112      !!          - prep = ridged and rafted ice + closing_gross
113      !!          - shift = move ice from one category to another
114      !!
115      !! ** Details
116      !!    step1: The net rate of closing is due to convergence and shear, based on Flato and Hibler (1995).
117      !!           The energy dissipation rate is equal to the net closing rate times the ice strength.
118      !!
119      !!    step3: The gross closing rate is equal to the first two terms (open
120      !!           water closing and thin ice ridging) without the third term
121      !!           (thick, newly ridged ice).
122      !!
123      !! References :   Flato, G. M., and W. D. Hibler III, 1995, JGR, 100, 18,611-18,626.
124      !!                Hibler, W. D. III, 1980, MWR, 108, 1943-1973, 1980.
125      !!                Rothrock, D. A., 1975: JGR, 80, 4514-4519.
126      !!                Thorndike et al., 1975, JGR, 80, 4501-4513.
127      !!                Bitz et al., JGR, 2001
128      !!                Amundrud and Melling, JGR 2005
129      !!                Babko et al., JGR 2002
130      !!
131      !!     This routine is based on CICE code and authors William H. Lipscomb,
132      !!     and Elizabeth C. Hunke, LANL are gratefully acknowledged
133      !!-------------------------------------------------------------------
134      INTEGER, INTENT(in) ::   kt     ! number of iteration
135      !!
136      INTEGER  ::   ji, jj, jk, jl             ! dummy loop index
137      INTEGER  ::   iter, iterate_ridging      ! local integer
138      INTEGER  ::   ipti                       ! local integer
139      REAL(wp) ::   zfac                       ! local scalar
140      INTEGER , DIMENSION(jpij) ::   iptidx        ! compute ridge/raft or not
141      REAL(wp), DIMENSION(jpij) ::   zdivu, zdelt  ! 1D divu_i & delta_i
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_2D( 1, 1, 1, 1 )
164         IF ( at_i(ji,jj) > epsi10 ) THEN
165            npti           = npti + 1
166            nptidx( npti ) = (jj - 1) * jpi + ji
167         ENDIF
168      END_2D
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_Dt_ice
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( sn_cfctl%l_prtctl )   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      ! In single precision there were floating point invalids due a sqrt of zhi which happens to have negative values
303      ! To solve that an extra check about the value of pv_i was added.
304      ! Although adding this condition is safe, the double definition (one for single other for double) has been kept to preserve the results of the sette test.
305#if defined key_single
306
307      WHERE( pa_i(1:npti,:) > epsi10 .and. pv_i(1:npti,:) > epsi10 )   ;   zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:)
308#else
309      WHERE( pa_i(1:npti,:) > epsi10 )   ;   zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:)
310#endif
311      ELSEWHERE                          ;   zhi(1:npti,:) = 0._wp
312      END WHERE
313
314      ! 1) Participation function (apartf): a(h) = b(h).g(h)
315      !-----------------------------------------------------------------
316      ! Compute the participation function = total area lost due to ridging/closing
317      ! This is analogous to
318      !   a(h) = b(h)g(h) as defined in Thorndike et al. (1975).
319      !   assuming b(h) = (2/Gstar) * (1 - G(h)/Gstar).
320      !
321      ! apartf = integrating b(h)g(h) between the category boundaries
322      ! apartf is always >= 0 and SUM(apartf(0:jpl))=1
323      !-----------------------------------------------------------------
324      !
325      ! Compute total area of ice plus open water.
326      ! This is in general not equal to one because of divergence during transport
327      zasum(1:npti) = pato_i(1:npti) + SUM( pa_i(1:npti,:), dim=2 )
328      !
329      WHERE( zasum(1:npti) > epsi10 )   ;   z1_asum(1:npti) = 1._wp / zasum(1:npti)
330      ELSEWHERE                         ;   z1_asum(1:npti) = 0._wp
331      END WHERE
332      !
333      ! Compute cumulative thickness distribution function
334      ! Compute the cumulative thickness distribution function zGsum,
335      ! where zGsum(n) is the fractional area in categories 0 to n.
336      ! initial value (in h = 0) = open water area
337      zGsum(1:npti,-1) = 0._wp
338      zGsum(1:npti,0 ) = pato_i(1:npti) * z1_asum(1:npti)
339      DO jl = 1, jpl
340         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)
341      END DO
342      !
343      IF( ln_partf_lin ) THEN          !--- Linear formulation (Thorndike et al., 1975)
344         DO jl = 0, jpl   
345            DO ji = 1, npti
346               IF    ( zGsum(ji,jl)   < rn_gstar ) THEN
347                  apartf(ji,jl) = z1_gstar * ( zGsum(ji,jl) - zGsum(ji,jl-1) ) * &
348                     &                       ( 2._wp - ( zGsum(ji,jl-1) + zGsum(ji,jl) ) * z1_gstar )
349               ELSEIF( zGsum(ji,jl-1) < rn_gstar ) THEN
350                  apartf(ji,jl) = z1_gstar * ( rn_gstar     - zGsum(ji,jl-1) ) *  &
351                     &                       ( 2._wp - ( zGsum(ji,jl-1) + rn_gstar        ) * z1_gstar )
352               ELSE
353                  apartf(ji,jl) = 0._wp
354               ENDIF
355            END DO
356         END DO
357         !
358      ELSEIF( ln_partf_exp ) THEN      !--- Exponential, more stable formulation (Lipscomb et al, 2007)
359         !                       
360         zfac = 1._wp / ( 1._wp - EXP(-z1_astar) )
361         DO jl = -1, jpl
362            DO ji = 1, npti
363               zGsum(ji,jl) = EXP( -zGsum(ji,jl) * z1_astar ) * zfac
364            END DO
365         END DO
366         DO jl = 0, jpl
367            DO ji = 1, npti
368               apartf(ji,jl) = zGsum(ji,jl-1) - zGsum(ji,jl)
369            END DO
370         END DO
371         !
372      ENDIF
373
374      !                                !--- Ridging and rafting participation concentrations
375      IF( ln_rafting .AND. ln_ridging ) THEN             !- ridging & rafting
376         DO jl = 1, jpl
377            DO ji = 1, npti
378               aridge(ji,jl) = ( 1._wp + TANH ( rn_craft * ( zhi(ji,jl) - rn_hraft ) ) ) * 0.5_wp * apartf(ji,jl)
379               araft (ji,jl) = apartf(ji,jl) - aridge(ji,jl)
380            END DO
381         END DO
382      ELSEIF( ln_ridging .AND. .NOT. ln_rafting ) THEN   !- ridging alone
383         DO jl = 1, jpl
384            DO ji = 1, npti
385               aridge(ji,jl) = apartf(ji,jl)
386               araft (ji,jl) = 0._wp
387            END DO
388         END DO
389      ELSEIF( ln_rafting .AND. .NOT. ln_ridging ) THEN   !- rafting alone   
390         DO jl = 1, jpl
391            DO ji = 1, npti
392               aridge(ji,jl) = 0._wp
393               araft (ji,jl) = apartf(ji,jl)
394            END DO
395         END DO
396      ELSE                                               !- no ridging & no rafting
397         DO jl = 1, jpl
398            DO ji = 1, npti
399               aridge(ji,jl) = 0._wp
400               araft (ji,jl) = 0._wp         
401            END DO
402         END DO
403      ENDIF
404
405      ! 2) Transfer function
406      !-----------------------------------------------------------------
407      ! Compute max and min ridged ice thickness for each ridging category.
408      ! Assume ridged ice is uniformly distributed between hrmin and hrmax.
409      !
410      ! This parameterization is a modified version of Hibler (1980).
411      ! The mean ridging thickness, zhmean, is proportional to hi^(0.5)
412      !  and for very thick ridging ice must be >= hrdg_hi_min*hi
413      !
414      ! The minimum ridging thickness, hrmin, is equal to 2*hi
415      !  (i.e., rafting) and for very thick ridging ice is
416      !  constrained by hrmin <= (zhmean + hi)/2.
417      !
418      ! The maximum ridging thickness, hrmax, is determined by zhmean and hrmin.
419      !
420      ! These modifications have the effect of reducing the ice strength
421      ! (relative to the Hibler formulation) when very thick ice is ridging.
422      !
423      ! zaksum = net area removed/ total area removed
424      ! where total area removed = area of ice that ridges
425      !         net area removed = total area removed - area of new ridges
426      !-----------------------------------------------------------------
427      zfac = 1._wp / hi_hrft
428      zaksum(1:npti) = apartf(1:npti,0)
429      !
430      DO jl = 1, jpl
431         DO ji = 1, npti
432            IF ( apartf(ji,jl) > 0._wp ) THEN
433               zhmean         = MAX( SQRT( rn_hstar * zhi(ji,jl) ), zhi(ji,jl) * hrdg_hi_min )
434               hrmin  (ji,jl) = MIN( 2._wp * zhi(ji,jl), 0.5_wp * ( zhmean + zhi(ji,jl) ) )
435               hrmax  (ji,jl) = 2._wp * zhmean - hrmin(ji,jl)
436               hraft  (ji,jl) = zhi(ji,jl) * zfac
437               hi_hrdg(ji,jl) = zhi(ji,jl) / MAX( zhmean, epsi20 )
438               !
439               ! Normalization factor : zaksum, ensures mass conservation
440               zaksum(ji) = zaksum(ji) + aridge(ji,jl) * ( 1._wp - hi_hrdg(ji,jl) )    &
441                  &                    + araft (ji,jl) * ( 1._wp - hi_hrft )
442            ELSE
443               hrmin  (ji,jl) = 0._wp 
444               hrmax  (ji,jl) = 0._wp 
445               hraft  (ji,jl) = 0._wp 
446               hi_hrdg(ji,jl) = 1._wp
447            ENDIF
448         END DO
449      END DO
450      !
451      ! 3) closing_gross
452      !-----------------
453      ! Based on the ITD of ridging and ridged ice, convert the net closing rate to a gross closing rate. 
454      ! NOTE: 0 < aksum <= 1
455      WHERE( zaksum(1:npti) > epsi10 )   ;   closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti)
456      ELSEWHERE                          ;   closing_gross(1:npti) = 0._wp
457      END WHERE
458     
459      ! correction to closing rate if excessive ice removal
460      !----------------------------------------------------
461      ! Reduce the closing rate if more than 100% of any ice category would be removed
462      ! Reduce the opening rate in proportion
463      DO jl = 1, jpl
464         DO ji = 1, npti
465            zfac = apartf(ji,jl) * closing_gross(ji) * rDt_ice
466            IF( zfac > pa_i(ji,jl) .AND. apartf(ji,jl) /= 0._wp ) THEN
467               closing_gross(ji) = pa_i(ji,jl) / apartf(ji,jl) * r1_Dt_ice
468            ENDIF
469         END DO
470      END DO     
471
472      ! 4) correction to opening if excessive open water removal
473      !---------------------------------------------------------
474      ! Reduce the closing rate if more than 100% of the open water would be removed
475      ! Reduce the opening rate in proportion
476      DO ji = 1, npti 
477         zfac = pato_i(ji) + ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * rDt_ice
478         IF( zfac < 0._wp ) THEN           ! would lead to negative ato_i
479            opning(ji) = apartf(ji,0) * closing_gross(ji) - pato_i(ji) * r1_Dt_ice 
480         ELSEIF( zfac > zasum(ji) ) THEN   ! would lead to ato_i > asum
481            opning(ji) = apartf(ji,0) * closing_gross(ji) + ( zasum(ji) - pato_i(ji) ) * r1_Dt_ice 
482         ENDIF
483      END DO
484      !
485   END SUBROUTINE rdgrft_prep
486
487
488   SUBROUTINE rdgrft_shift
489      !!-------------------------------------------------------------------
490      !!                ***  ROUTINE rdgrft_shift ***
491      !!
492      !! ** Purpose :   shift ridging ice among thickness categories of ice thickness
493      !!
494      !! ** Method  :   Remove area, volume, and energy from each ridging category
495      !!                and add to thicker ice categories.
496      !!-------------------------------------------------------------------
497      !
498      INTEGER  ::   ji, jj, jl, jl1, jl2, jk   ! dummy loop indices
499      REAL(wp) ::   hL, hR, farea              ! left and right limits of integration and new area going to jl2
500      REAL(wp) ::   vsw                        ! vol of water trapped into ridges
501      REAL(wp) ::   afrdg, afrft               ! fraction of category area ridged/rafted
502      REAL(wp)                  ::   airdg1, oirdg1, aprdg1, virdg1, sirdg1
503      REAL(wp)                  ::   airft1, oirft1, aprft1
504      REAL(wp), DIMENSION(jpij) ::   airdg2, oirdg2, aprdg2, virdg2, sirdg2, vsrdg, vprdg, vlrdg  ! area etc of new ridges
505      REAL(wp), DIMENSION(jpij) ::   airft2, oirft2, aprft2, virft , sirft , vsrft, vprft, vlrft  ! area etc of rafted ice
506      !
507      REAL(wp), DIMENSION(jpij) ::   ersw             ! enth of water trapped into ridges
508      REAL(wp), DIMENSION(jpij) ::   zswitch, fvol    ! new ridge volume going to jl2
509      REAL(wp), DIMENSION(jpij) ::   z1_ai            ! 1 / a
510      REAL(wp), DIMENSION(jpij) ::   zvti             ! sum(v_i)
511      !
512      REAL(wp), DIMENSION(jpij,nlay_s) ::   esrft     ! snow energy of rafting ice
513      REAL(wp), DIMENSION(jpij,nlay_i) ::   eirft     ! ice  energy of rafting ice
514      REAL(wp), DIMENSION(jpij,nlay_s) ::   esrdg     ! enth*volume of new ridges     
515      REAL(wp), DIMENSION(jpij,nlay_i) ::   eirdg     ! enth*volume of new ridges
516      !
517      INTEGER , DIMENSION(jpij) ::   itest_rdg, itest_rft   ! test for conservation
518      !!-------------------------------------------------------------------
519      !
520      zvti(1:npti) = SUM( v_i_2d(1:npti,:), dim=2 )   ! total ice volume
521      !
522      ! 1) Change in open water area due to closing and opening
523      !--------------------------------------------------------
524      DO ji = 1, npti
525         ato_i_1d(ji) = MAX( 0._wp, ato_i_1d(ji) + ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * rDt_ice )
526      END DO
527     
528      ! 2) compute categories in which ice is removed (jl1)
529      !----------------------------------------------------
530      DO jl1 = 1, jpl
531
532         IF( nn_icesal /= 2 )  THEN     
533            CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d(1:npti), s_i(:,:,jl1) )
534         ENDIF
535
536         DO ji = 1, npti
537
538            IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN   ! only if ice is ridging
539
540               IF( a_i_2d(ji,jl1) > epsi10 ) THEN   ;   z1_ai(ji) = 1._wp / a_i_2d(ji,jl1)
541               ELSE                                 ;   z1_ai(ji) = 0._wp
542               ENDIF
543               
544               ! area of ridging / rafting ice (airdg1) and of new ridge (airdg2)
545               airdg1 = aridge(ji,jl1) * closing_gross(ji) * rDt_ice
546               airft1 = araft (ji,jl1) * closing_gross(ji) * rDt_ice
547
548               airdg2(ji) = airdg1 * hi_hrdg(ji,jl1)
549               airft2(ji) = airft1 * hi_hrft
550
551               ! ridging /rafting fractions
552               afrdg = airdg1 * z1_ai(ji)
553               afrft = airft1 * z1_ai(ji)
554
555               ! volume and enthalpy (J/m2, >0) of seawater trapped into ridges
556               IF    ( zvti(ji) <= 10. ) THEN ; vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg                                           ! v <= 10m then porosity = rn_porordg
557               ELSEIF( zvti(ji) >= 20. ) THEN ; vsw = 0._wp                                                                         ! v >= 20m then porosity = 0
558               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
559               ENDIF
560               ersw(ji) = -rhoi * vsw * rcp * sst_1d(ji)   ! clem: if sst>0, then ersw <0 (is that possible?)
561
562               ! volume etc of ridging / rafting ice and new ridges (vi, vs, sm, oi, es, ei)
563               virdg1     = v_i_2d (ji,jl1)   * afrdg
564               virdg2(ji) = v_i_2d (ji,jl1)   * afrdg + vsw
565               vsrdg(ji)  = v_s_2d (ji,jl1)   * afrdg
566               sirdg1     = sv_i_2d(ji,jl1)   * afrdg
567               sirdg2(ji) = sv_i_2d(ji,jl1)   * afrdg + vsw * sss_1d(ji)
568               oirdg1     = oa_i_2d(ji,jl1)   * afrdg
569               oirdg2(ji) = oa_i_2d(ji,jl1)   * afrdg * hi_hrdg(ji,jl1) 
570
571               virft(ji)  = v_i_2d (ji,jl1)   * afrft
572               vsrft(ji)  = v_s_2d (ji,jl1)   * afrft
573               sirft(ji)  = sv_i_2d(ji,jl1)   * afrft 
574               oirft1     = oa_i_2d(ji,jl1)   * afrft 
575               oirft2(ji) = oa_i_2d(ji,jl1)   * afrft * hi_hrft 
576
577               IF ( ln_pnd_LEV ) THEN
578                  aprdg1     = a_ip_2d(ji,jl1) * afrdg
579                  aprdg2(ji) = a_ip_2d(ji,jl1) * afrdg * hi_hrdg(ji,jl1)
580                  vprdg (ji) = v_ip_2d(ji,jl1) * afrdg
581                  aprft1     = a_ip_2d(ji,jl1) * afrft
582                  aprft2(ji) = a_ip_2d(ji,jl1) * afrft * hi_hrft
583                  vprft (ji) = v_ip_2d(ji,jl1) * afrft
584                  IF ( ln_pnd_lids ) THEN
585                     vlrdg (ji) = v_il_2d(ji,jl1) * afrdg
586                     vlrft (ji) = v_il_2d(ji,jl1) * afrft
587                  ENDIF
588               ENDIF
589
590               ! Ice-ocean exchanges associated with ice porosity
591               wfx_dyn_1d(ji) = wfx_dyn_1d(ji) - vsw * rhoi * r1_Dt_ice   ! increase in ice volume due to seawater frozen in voids
592               sfx_dyn_1d(ji) = sfx_dyn_1d(ji) - vsw * sss_1d(ji) * rhoi * r1_Dt_ice
593               hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ersw(ji) * r1_Dt_ice          ! > 0 [W.m-2]
594
595               ! Put the snow lost by ridging into the ocean
596               !  Note that esrdg > 0; the ocean must cool to melt snow. If the ocean temp = Tf already, new ice must grow.
597               wfx_snw_dyn_1d(ji) = wfx_snw_dyn_1d(ji) + ( rhos * vsrdg(ji) * ( 1._wp - rn_fsnwrdg )   &   ! fresh water source for ocean
598                  &                                      + rhos * vsrft(ji) * ( 1._wp - rn_fsnwrft ) ) * r1_Dt_ice
599
600               ! virtual salt flux to keep salinity constant
601               IF( nn_icesal /= 2 )  THEN
602                  sirdg2(ji)     = sirdg2(ji)     - vsw * ( sss_1d(ji) - s_i_1d(ji) )       ! ridge salinity = s_i
603                  sfx_bri_1d(ji) = sfx_bri_1d(ji) + sss_1d(ji) * vsw * rhoi * r1_Dt_ice  &  ! put back sss_m into the ocean
604                     &                            - s_i_1d(ji) * vsw * rhoi * r1_Dt_ice     ! and get  s_i  from the ocean
605               ENDIF
606
607               ! Remove area, volume of new ridge to each category jl1
608               !------------------------------------------------------
609               a_i_2d (ji,jl1) = a_i_2d (ji,jl1) - airdg1    - airft1
610               v_i_2d (ji,jl1) = v_i_2d (ji,jl1) - virdg1    - virft(ji)
611               v_s_2d (ji,jl1) = v_s_2d (ji,jl1) - vsrdg(ji) - vsrft(ji)
612               sv_i_2d(ji,jl1) = sv_i_2d(ji,jl1) - sirdg1    - sirft(ji)
613               oa_i_2d(ji,jl1) = oa_i_2d(ji,jl1) - oirdg1    - oirft1
614               IF ( ln_pnd_LEV ) THEN
615                  a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - aprdg1    - aprft1
616                  v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - vprdg(ji) - vprft(ji)
617                  IF ( ln_pnd_lids ) THEN
618                     v_il_2d(ji,jl1) = v_il_2d(ji,jl1) - vlrdg(ji) - vlrft(ji)
619                  ENDIF
620               ENDIF
621            ENDIF
622
623         END DO ! ji
624
625         ! special loop for e_s because of layers jk
626         DO jk = 1, nlay_s
627            DO ji = 1, npti
628               IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
629                  ! Compute ridging /rafting fractions
630                  afrdg = aridge(ji,jl1) * closing_gross(ji) * rDt_ice * z1_ai(ji)
631                  afrft = araft (ji,jl1) * closing_gross(ji) * rDt_ice * z1_ai(ji)
632                  ! Compute ridging /rafting ice and new ridges for es
633                  esrdg(ji,jk) = ze_s_2d (ji,jk,jl1) * afrdg
634                  esrft(ji,jk) = ze_s_2d (ji,jk,jl1) * afrft
635                  ! Put the snow lost by ridging into the ocean
636                  hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ( - esrdg(ji,jk) * ( 1._wp - rn_fsnwrdg )   &                 ! heat sink for ocean (<0, W.m-2)
637                     &                                - esrft(ji,jk) * ( 1._wp - rn_fsnwrft ) ) * r1_Dt_ice
638                  !
639                  ! Remove energy of new ridge to each category jl1
640                  !-------------------------------------------------
641                  ze_s_2d(ji,jk,jl1) = ze_s_2d(ji,jk,jl1) * ( 1._wp - afrdg - afrft ) 
642               ENDIF
643            END DO
644         END DO
645                 
646         ! special loop for e_i because of layers jk
647         DO jk = 1, nlay_i
648            DO ji = 1, npti
649               IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
650                  ! Compute ridging /rafting fractions
651                  afrdg = aridge(ji,jl1) * closing_gross(ji) * rDt_ice * z1_ai(ji)
652                  afrft = araft (ji,jl1) * closing_gross(ji) * rDt_ice * z1_ai(ji)
653                  ! Compute ridging ice and new ridges for ei
654                  eirdg(ji,jk) = ze_i_2d (ji,jk,jl1) * afrdg + ersw(ji) * r1_nlay_i
655                  eirft(ji,jk) = ze_i_2d (ji,jk,jl1) * afrft
656                  !
657                  ! Remove energy of new ridge to each category jl1
658                  !-------------------------------------------------
659                  ze_i_2d(ji,jk,jl1) = ze_i_2d(ji,jk,jl1) * ( 1._wp - afrdg - afrft ) 
660               ENDIF
661            END DO
662         END DO
663         
664         ! 3) compute categories in which ice is added (jl2)
665         !--------------------------------------------------
666         itest_rdg(1:npti) = 0
667         itest_rft(1:npti) = 0
668         DO jl2  = 1, jpl 
669            !
670            DO ji = 1, npti
671
672               IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
673
674                  ! Compute the fraction of ridged ice area and volume going to thickness category jl2
675                  IF( hrmin(ji,jl1) <= hi_max(jl2) .AND. hrmax(ji,jl1) > hi_max(jl2-1) ) THEN
676                     hL = MAX( hrmin(ji,jl1), hi_max(jl2-1) )
677                     hR = MIN( hrmax(ji,jl1), hi_max(jl2)   )
678                     farea    = ( hR      - hL      ) / ( hrmax(ji,jl1)                 - hrmin(ji,jl1)                 )
679                     fvol(ji) = ( hR * hR - hL * hL ) / ( hrmax(ji,jl1) * hrmax(ji,jl1) - hrmin(ji,jl1) * hrmin(ji,jl1) )
680                     !
681                     itest_rdg(ji) = 1   ! test for conservation
682                  ELSE
683                     farea    = 0._wp 
684                     fvol(ji) = 0._wp                 
685                  ENDIF
686
687                  ! Compute the fraction of rafted ice area and volume going to thickness category jl2
688                  IF( hraft(ji,jl1) <= hi_max(jl2) .AND. hraft(ji,jl1) >  hi_max(jl2-1) ) THEN
689                     zswitch(ji) = 1._wp
690                     !
691                     itest_rft(ji) = 1   ! test for conservation
692                  ELSE
693                     zswitch(ji) = 0._wp
694                  ENDIF
695                  !
696                  ! Patch to ensure perfect conservation if ice thickness goes mad
697                  ! Sometimes thickness is larger than hi_max(jpl) because of advection scheme (for very small areas)
698                  ! Then ice volume is removed from one category but the ridging/rafting scheme
699                  ! does not know where to move it, leading to a conservation issue. 
700                  IF( itest_rdg(ji) == 0 .AND. jl2 == jpl ) THEN   ;   farea = 1._wp   ;   fvol(ji) = 1._wp   ;   ENDIF
701                  IF( itest_rft(ji) == 0 .AND. jl2 == jpl )      zswitch(ji) = 1._wp
702                  !
703                  ! Add area, volume of new ridge to category jl2
704                  !----------------------------------------------
705                  a_i_2d (ji,jl2) = a_i_2d (ji,jl2) + ( airdg2(ji) * farea    + airft2(ji) * zswitch(ji) )
706                  oa_i_2d(ji,jl2) = oa_i_2d(ji,jl2) + ( oirdg2(ji) * farea    + oirft2(ji) * zswitch(ji) )
707                  v_i_2d (ji,jl2) = v_i_2d (ji,jl2) + ( virdg2(ji) * fvol(ji) + virft (ji) * zswitch(ji) )
708                  sv_i_2d(ji,jl2) = sv_i_2d(ji,jl2) + ( sirdg2(ji) * fvol(ji) + sirft (ji) * zswitch(ji) )
709                  v_s_2d (ji,jl2) = v_s_2d (ji,jl2) + ( vsrdg (ji) * rn_fsnwrdg * fvol(ji)  +  &
710                     &                                  vsrft (ji) * rn_fsnwrft * zswitch(ji) )
711                  IF ( ln_pnd_LEV ) THEN
712                     v_ip_2d (ji,jl2) = v_ip_2d(ji,jl2) + (   vprdg (ji) * rn_fpndrdg * fvol   (ji)   &
713                        &                                   + vprft (ji) * rn_fpndrft * zswitch(ji)   )
714                     a_ip_2d (ji,jl2) = a_ip_2d(ji,jl2) + (   aprdg2(ji) * rn_fpndrdg * farea         & 
715                        &                                   + aprft2(ji) * rn_fpndrft * zswitch(ji)   )
716                     IF ( ln_pnd_lids ) THEN
717                        v_il_2d (ji,jl2) = v_il_2d(ji,jl2) + (   vlrdg(ji) * rn_fpndrdg * fvol   (ji) &
718                           &                                   + vlrft(ji) * rn_fpndrft * zswitch(ji) )
719                     ENDIF
720                  ENDIF
721                 
722               ENDIF
723
724            END DO
725            ! Add snow energy of new ridge to category jl2
726            !---------------------------------------------
727            DO jk = 1, nlay_s
728               DO ji = 1, npti
729                  IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp )   &
730                     &   ze_s_2d(ji,jk,jl2) = ze_s_2d(ji,jk,jl2) + ( esrdg(ji,jk) * rn_fsnwrdg * fvol(ji)  +  &
731                     &                                               esrft(ji,jk) * rn_fsnwrft * zswitch(ji) )
732               END DO
733            END DO
734            ! Add ice energy of new ridge to category jl2
735            !--------------------------------------------
736            DO jk = 1, nlay_i
737               DO ji = 1, npti
738                  IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp )   &
739                     &   ze_i_2d(ji,jk,jl2) = ze_i_2d(ji,jk,jl2) + eirdg(ji,jk) * fvol(ji) + eirft(ji,jk) * zswitch(ji)                 
740               END DO
741            END DO
742            !
743         END DO ! jl2
744         !
745      END DO ! jl1
746      !
747      ! roundoff errors
748      !----------------
749      ! In case ridging/rafting lead to very small negative values (sometimes it happens)
750      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 )
751      !
752   END SUBROUTINE rdgrft_shift
753
754
755   SUBROUTINE ice_strength
756      !!----------------------------------------------------------------------
757      !!                ***  ROUTINE ice_strength ***
758      !!
759      !! ** Purpose :   computes ice strength used in dynamics routines of ice thickness
760      !!
761      !! ** Method  :   Compute the strength of the ice pack, defined as the energy (J m-2)
762      !!              dissipated per unit area removed from the ice pack under compression,
763      !!              and assumed proportional to the change in potential energy caused
764      !!              by ridging. Note that only Hibler's formulation is stable and that
765      !!              ice strength has to be smoothed
766      !!----------------------------------------------------------------------
767      INTEGER             ::   ji, jj, jl  ! dummy loop indices
768      INTEGER             ::   ismooth     ! smoothing the resistance to deformation
769      INTEGER             ::   itframe     ! number of time steps for the P smoothing
770      REAL(wp)            ::   zp, z1_3    ! local scalars
771      REAL(wp), DIMENSION(jpi,jpj) ::   zworka           ! temporary array used here
772      REAL(wp), DIMENSION(jpi,jpj) ::   zstrp1, zstrp2   ! strength at previous time steps
773      !!----------------------------------------------------------------------
774      !                              !--------------------------------------------------!
775      IF( ln_str_H79 ) THEN          ! Ice strength => Hibler (1979) method             !
776      !                              !--------------------------------------------------!
777         strength(:,:) = rn_pstar * SUM( v_i(:,:,:), dim=3 ) * EXP( -rn_crhg * ( 1._wp - SUM( a_i(:,:,:), dim=3 ) ) )
778         ismooth = 1
779         !                           !--------------------------------------------------!
780      ELSE                           ! Zero strength                                    !
781         !                           !--------------------------------------------------!
782         strength(:,:) = 0._wp
783         ismooth = 0
784      ENDIF
785      !                              !--------------------------------------------------!
786      SELECT CASE( ismooth )         ! Smoothing ice strength                           !
787      !                              !--------------------------------------------------!
788      CASE( 1 )               !--- Spatial smoothing
789         DO_2D( 0, 0, 0, 0 )
790            IF ( SUM( a_i(ji,jj,:) ) > 0._wp ) THEN
791               zworka(ji,jj) = ( 4.0 * strength(ji,jj)              &
792                  &                  + strength(ji-1,jj) * tmask(ji-1,jj,1) + strength(ji+1,jj) * tmask(ji+1,jj,1) & 
793                  &                  + strength(ji,jj-1) * tmask(ji,jj-1,1) + strength(ji,jj+1) * tmask(ji,jj+1,1) &
794                  &            ) / ( 4.0 + tmask(ji-1,jj,1) + tmask(ji+1,jj,1) + tmask(ji,jj-1,1) + tmask(ji,jj+1,1) )
795            ELSE
796               zworka(ji,jj) = 0._wp
797            ENDIF
798         END_2D
799         
800         DO_2D( 0, 0, 0, 0 )
801            strength(ji,jj) = zworka(ji,jj)
802         END_2D
803         CALL lbc_lnk( 'icedyn_rdgrft', strength, 'T', 1.0_wp )
804         !
805      CASE( 2 )               !--- Temporal smoothing
806         IF ( kt_ice == nit000 ) THEN
807            zstrp1(:,:) = 0._wp
808            zstrp2(:,:) = 0._wp
809         ENDIF
810         !
811         DO_2D( 0, 0, 0, 0 )
812            IF ( SUM( a_i(ji,jj,:) ) > 0._wp ) THEN
813               itframe = 1 ! number of time steps for the running mean
814               IF ( zstrp1(ji,jj) > 0._wp ) itframe = itframe + 1
815               IF ( zstrp2(ji,jj) > 0._wp ) itframe = itframe + 1
816               zp = ( strength(ji,jj) + zstrp1(ji,jj) + zstrp2(ji,jj) ) / itframe
817               zstrp2  (ji,jj) = zstrp1  (ji,jj)
818               zstrp1  (ji,jj) = strength(ji,jj)
819               strength(ji,jj) = zp
820            ENDIF
821         END_2D
822         CALL lbc_lnk( 'icedyn_rdgrft', strength, 'T', 1.0_wp )
823         !
824      END SELECT
825      !
826   END SUBROUTINE ice_strength
827
828   
829   SUBROUTINE ice_dyn_1d2d( kn )
830      !!-----------------------------------------------------------------------
831      !!                   ***  ROUTINE ice_dyn_1d2d ***
832      !!                 
833      !! ** Purpose :   move arrays from 1d to 2d and the reverse
834      !!-----------------------------------------------------------------------
835      INTEGER, INTENT(in) ::   kn   ! 1= from 2D to 1D   ;   2= from 1D to 2D
836      !
837      INTEGER ::   jl, jk   ! dummy loop indices
838      !!-----------------------------------------------------------------------
839      !
840      SELECT CASE( kn )
841      !                    !---------------------!
842      CASE( 1 )            !==  from 2D to 1D  ==!
843         !                 !---------------------!
844         ! fields used but not modified
845         CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d(1:npti), sss_m(:,:) )
846         CALL tab_2d_1d( npti, nptidx(1:npti), sst_1d(1:npti), sst_m(:,:) )
847         ! the following fields are modified in this routine
848         !!CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti), ato_i(:,:) )
849         !!CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d(1:npti,1:jpl), a_i(:,:,:) )
850         !!CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d  (1:npti,1:jpl), v_i  (:,:,:) )
851         CALL tab_3d_2d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s (:,:,:) )
852         CALL tab_3d_2d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
853         CALL tab_3d_2d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i(:,:,:) )
854         CALL tab_3d_2d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) )
855         CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) )
856         CALL tab_3d_2d( npti, nptidx(1:npti), v_il_2d(1:npti,1:jpl), v_il(:,:,:) )
857         DO jl = 1, jpl
858            DO jk = 1, nlay_s
859               CALL tab_2d_1d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) )
860            END DO
861            DO jk = 1, nlay_i
862               CALL tab_2d_1d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
863            END DO
864         END DO
865         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_dyn_1d    (1:npti), sfx_dyn    (:,:) )
866         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_bri_1d    (1:npti), sfx_bri    (:,:) )
867         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_dyn_1d    (1:npti), wfx_dyn    (:,:) )
868         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_dyn_1d    (1:npti), hfx_dyn    (:,:) )
869         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) )
870         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d    (1:npti), wfx_pnd    (:,:) )
871         !
872         !                 !---------------------!
873      CASE( 2 )            !==  from 1D to 2D  ==!
874         !                 !---------------------!
875         CALL tab_1d_2d( npti, nptidx(1:npti), ato_i_1d(1:npti), ato_i(:,:) )
876         CALL tab_2d_3d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i (:,:,:) )
877         CALL tab_2d_3d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) )
878         CALL tab_2d_3d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s (:,:,:) )
879         CALL tab_2d_3d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
880         CALL tab_2d_3d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i(:,:,:) )
881         CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) )
882         CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) )
883         CALL tab_2d_3d( npti, nptidx(1:npti), v_il_2d(1:npti,1:jpl), v_il(:,:,:) )
884         DO jl = 1, jpl
885            DO jk = 1, nlay_s
886               CALL tab_1d_2d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) )
887            END DO
888            DO jk = 1, nlay_i
889               CALL tab_1d_2d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
890            END DO
891         END DO
892         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_dyn_1d    (1:npti), sfx_dyn    (:,:) )
893         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_bri_1d    (1:npti), sfx_bri    (:,:) )
894         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_dyn_1d    (1:npti), wfx_dyn    (:,:) )
895         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_dyn_1d    (1:npti), hfx_dyn    (:,:) )
896         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) )
897         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d    (1:npti), wfx_pnd    (:,:) )
898         !
899      END SELECT
900      !
901   END SUBROUTINE ice_dyn_1d2d
902   
903
904   SUBROUTINE ice_dyn_rdgrft_init
905      !!-------------------------------------------------------------------
906      !!                  ***  ROUTINE ice_dyn_rdgrft_init ***
907      !!
908      !! ** Purpose :   Physical constants and parameters linked
909      !!                to the mechanical ice redistribution
910      !!
911      !! ** Method  :   Read the namdyn_rdgrft namelist
912      !!                and check the parameters values
913      !!                called at the first timestep (nit000)
914      !!
915      !! ** input   :   Namelist namdyn_rdgrft
916      !!-------------------------------------------------------------------
917      INTEGER :: ios                 ! Local integer output status for namelist read
918      !!
919      NAMELIST/namdyn_rdgrft/ ln_str_H79, rn_pstar, rn_crhg, &
920         &                    rn_csrdg  ,                    &
921         &                    ln_partf_lin, rn_gstar,        &
922         &                    ln_partf_exp, rn_astar,        & 
923         &                    ln_ridging, rn_hstar, rn_porordg, rn_fsnwrdg, rn_fpndrdg,  & 
924         &                    ln_rafting, rn_hraft, rn_craft  , rn_fsnwrft, rn_fpndrft
925      !!-------------------------------------------------------------------
926      !
927      READ  ( numnam_ice_ref, namdyn_rdgrft, IOSTAT = ios, ERR = 901)
928901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdyn_rdgrft in reference namelist' )
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.