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 @ 15344

Last change on this file since 15344 was 15344, checked in by emmafiedler, 3 years ago

Include Rothrock ice strength formulation

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