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

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

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