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/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE – NEMO

source: NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icedyn_rdgrft.F90 @ 11501

Last change on this file since 11501 was 11501, checked in by clem, 5 years ago

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