New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
icedyn_rdgrft.F90 in NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/ICE – NEMO

source: NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/ICE/icedyn_rdgrft.F90 @ 14015

Last change on this file since 14015 was 14015, checked in by rlod, 4 years ago

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