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/UKMO/NEMO_4.0.4_ice_strength/src/ICE – NEMO

source: NEMO/branches/UKMO/NEMO_4.0.4_ice_strength/src/ICE/icedyn_rdgrft.F90 @ 15542

Last change on this file since 15542 was 15542, checked in by emmafiedler, 8 months ago

Exponential ridge ice redistribution

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