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

source: NEMO/branches/2021/ticket2632_r14588_theta_sbcblk/src/ICE/icedyn_rdgrft.F90 @ 15548

Last change on this file since 15548 was 15548, checked in by gsamson, 3 years ago

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