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_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/ICE – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/ICE/icedyn_rdgrft.F90 @ 10998

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

remove nn_virtual_itd options and replace them by ln_virtual_itd
  • Property svn:keywords set to Id
File size: 50.5 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      ! clem: The redistribution of ice between categories can lead to small negative values (as for the remapping in ice_itd_rem)
145      !       likely due to truncation error ( i.e. 1. - 1. /= 0 )
146      !       I do not think it should be a concern since small areas and volumes are erased (in ice_var_zapsmall.F90)
147     
148      ! controls
149      IF( ln_timing    )   CALL timing_start('icedyn_rdgrft')                                                             ! timing
150      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icedyn_rdgrft', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
151
152      IF( kt == nit000 ) THEN
153         IF(lwp) WRITE(numout,*)
154         IF(lwp) WRITE(numout,*)'ice_dyn_rdgrft: ice ridging and rafting'
155         IF(lwp) WRITE(numout,*)'~~~~~~~~~~~~~~'
156      ENDIF     
157
158      CALL ice_var_zapsmall   ! Zero out categories with very small areas
159
160      !--------------------------------
161      ! 0) Identify grid cells with ice
162      !--------------------------------
163      npti = 0   ;   nptidx(:) = 0
164      ipti = 0   ;   iptidx(:) = 0
165      DO jj = 1, jpj
166         DO ji = 1, jpi
167            IF ( at_i(ji,jj) > 0._wp ) THEN
168               npti           = npti + 1
169               nptidx( npti ) = (jj - 1) * jpi + ji
170            ENDIF
171         END DO
172      END DO
173     
174      !--------------------------------------------------------
175      ! 1) Dynamical inputs (closing rate, divergence, opening)
176      !--------------------------------------------------------
177      IF( npti > 0 ) THEN
178       
179         ! just needed here
180         CALL tab_2d_1d( npti, nptidx(1:npti), zdivu(1:npti), divu_i(:,:) )
181         CALL tab_2d_1d( npti, nptidx(1:npti), zdelt(1:npti), delta_i(:,:) )
182         ! needed here and in the iteration loop
183         CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d  (1:npti,1:jpl), a_i(:,:,:) )
184         CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d  (1:npti,1:jpl), v_i(:,:,:) )
185         CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti)      , ato_i(:,:) )
186
187         DO ji = 1, npti
188            ! closing_net = rate at which open water area is removed + ice area removed by ridging
189            !                                                        - ice area added in new ridges
190            closing_net(ji) = rn_csrdg * 0.5_wp * ( zdelt(ji) - ABS( zdivu(ji) ) ) - MIN( zdivu(ji), 0._wp )
191            !
192            ! divergence given by the advection scheme
193            !   (which may not be equal to divu as computed from the velocity field)
194            IF    ( ln_adv_Pra ) THEN
195               zdivu_adv(ji) = ( 1._wp - ato_i_1d(ji) - SUM( a_i_2d(ji,:) ) ) * r1_rdtice
196            ELSEIF( ln_adv_UMx ) THEN
197               zdivu_adv(ji) = zdivu(ji)
198            ENDIF
199            !
200            IF( zdivu_adv(ji) < 0._wp )   closing_net(ji) = MAX( closing_net(ji), -zdivu_adv(ji) )   ! make sure the closing rate is large enough
201            !                                                                                        ! to give asum = 1.0 after ridging
202            ! Opening rate (non-negative) that will give asum = 1.0 after ridging.
203            opning(ji) = closing_net(ji) + zdivu_adv(ji)
204         END DO
205         !
206         !------------------------------------
207         ! 2) Identify grid cells with ridging
208         !------------------------------------
209         CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d, closing_net )
210
211         DO ji = 1, npti
212            IF( SUM( apartf(ji,1:jpl) ) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
213               ipti = ipti + 1
214               iptidx     (ipti)   = nptidx     (ji)
215               ! adjust to new indices
216               a_i_2d     (ipti,:) = a_i_2d     (ji,:)
217               v_i_2d     (ipti,:) = v_i_2d     (ji,:)
218               ato_i_1d   (ipti)   = ato_i_1d   (ji)
219               closing_net(ipti)   = closing_net(ji)
220               zdivu_adv  (ipti)   = zdivu_adv  (ji)
221               opning     (ipti)   = opning     (ji)
222            ENDIF
223         END DO
224
225      ENDIF
226
227      ! grid cells with ridging
228      nptidx(:) = iptidx(:)
229      npti      = ipti
230
231      !-----------------
232      ! 3) Start ridging
233      !-----------------
234      IF( npti > 0 ) THEN
235         
236         CALL ice_dyn_1d2d( 1 )            ! --- Move to 1D arrays --- !
237
238         iter            = 1
239         iterate_ridging = 1     
240         !                                                        !----------------------!
241         DO WHILE( iterate_ridging > 0 .AND. iter < jp_itermax )  !  ridging iterations  !
242            !                                                     !----------------------!
243            ! Calculate participation function (apartf)
244            !       and transfer      function
245            !       and closing_gross (+correction on opening)
246            CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d, closing_net )
247
248            ! Redistribute area, volume, and energy between categories
249            CALL rdgrft_shift
250
251            ! Do we keep on iterating?
252            !-------------------------
253            ! Check whether a_i + ato_i = 0
254            ! If not, because the closing and opening rates were reduced above, ridge again with new rates
255            iterate_ridging = 0
256            DO ji = 1, npti
257               zfac = 1._wp - ( ato_i_1d(ji) + SUM( a_i_2d(ji,:) ) )
258               IF( ABS( zfac ) < epsi10 ) THEN
259                  closing_net(ji) = 0._wp
260                  opning     (ji) = 0._wp
261                  ato_i_1d   (ji) = MAX( 0._wp, 1._wp - SUM( a_i_2d(ji,:) ) )
262               ELSE
263                  iterate_ridging  = 1
264                  zdivu_adv  (ji) = zfac * r1_rdtice
265                  closing_net(ji) = MAX( 0._wp, -zdivu_adv(ji) )
266                  opning     (ji) = MAX( 0._wp,  zdivu_adv(ji) )
267               ENDIF
268            END DO
269            !
270            iter = iter + 1
271            IF( iter  >  jp_itermax )    CALL ctl_stop( 'STOP',  'icedyn_rdgrft: non-converging ridging scheme'  )
272            !
273         END DO
274
275         CALL ice_dyn_1d2d( 2 )            ! --- Move to 2D arrays --- !
276
277      ENDIF
278   
279      CALL ice_var_agg( 1 ) 
280
281      ! controls
282      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icedyn_rdgrft', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
283      IF( ln_ctl       )   CALL ice_prt3D   ('icedyn_rdgrft')                                                             ! prints
284      IF( ln_timing    )   CALL timing_stop ('icedyn_rdgrft')                                                             ! timing
285      !
286   END SUBROUTINE ice_dyn_rdgrft
287
288
289   SUBROUTINE rdgrft_prep( pa_i, pv_i, pato_i, pclosing_net )
290      !!-------------------------------------------------------------------
291      !!                ***  ROUTINE rdgrft_prep ***
292      !!
293      !! ** Purpose :   preparation for ridging calculations
294      !!
295      !! ** Method  :   Compute the thickness distribution of the ice and open water
296      !!                participating in ridging and of the resulting ridges.
297      !!-------------------------------------------------------------------
298      REAL(wp), DIMENSION(:)  , INTENT(in) ::   pato_i, pclosing_net 
299      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pa_i, pv_i 
300      !!
301      INTEGER  ::   ji, jl                     ! dummy loop indices
302      REAL(wp) ::   z1_gstar, z1_astar, zhmean, zfac   ! local scalar
303      REAL(wp), DIMENSION(jpij)        ::   zasum, z1_asum, zaksum   ! sum of a_i+ato_i and reverse
304      REAL(wp), DIMENSION(jpij,jpl)    ::   zhi                      ! ice thickness
305      REAL(wp), DIMENSION(jpij,-1:jpl) ::   zGsum                    ! zGsum(n) = sum of areas in categories 0 to n
306      !--------------------------------------------------------------------
307
308      z1_gstar = 1._wp / rn_gstar
309      z1_astar = 1._wp / rn_astar
310
311      !                       ! Ice thickness needed for rafting
312      WHERE( pa_i(1:npti,:) > 0._wp )   ;   zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:)
313      ELSEWHERE                         ;   zhi(1:npti,:) = 0._wp
314      END WHERE
315
316      ! 1) Participation function (apartf): a(h) = b(h).g(h)
317      !-----------------------------------------------------------------
318      ! Compute the participation function = total area lost due to ridging/closing
319      ! This is analogous to
320      !   a(h) = b(h)g(h) as defined in Thorndike et al. (1975).
321      !   assuming b(h) = (2/Gstar) * (1 - G(h)/Gstar).
322      !
323      ! apartf = integrating b(h)g(h) between the category boundaries
324      ! apartf is always >= 0 and SUM(apartf(0:jpl))=1
325      !-----------------------------------------------------------------
326      !
327      ! Compute total area of ice plus open water.
328      ! This is in general not equal to one because of divergence during transport
329      zasum(1:npti) = pato_i(1:npti) + SUM( pa_i(1:npti,:), dim=2 )
330      !
331      WHERE( zasum(1:npti) > 0._wp )   ;   z1_asum(1:npti) = 1._wp / zasum(1:npti)
332      ELSEWHERE                        ;   z1_asum(1:npti) = 0._wp
333      END WHERE
334      !
335      ! Compute cumulative thickness distribution function
336      ! Compute the cumulative thickness distribution function zGsum,
337      ! where zGsum(n) is the fractional area in categories 0 to n.
338      ! initial value (in h = 0) = open water area
339      zGsum(1:npti,-1) = 0._wp
340      zGsum(1:npti,0 ) = pato_i(1:npti) * z1_asum(1:npti)
341      DO jl = 1, jpl
342         zGsum(1:npti,jl) = ( pato_i(1:npti) + SUM( pa_i(1:npti,1:jl), dim=2 ) ) * z1_asum(1:npti)  ! sum(1:jl) is ok (and not jpl)
343      END DO
344      !
345      IF( ln_partf_lin ) THEN          !--- Linear formulation (Thorndike et al., 1975)
346         DO jl = 0, jpl   
347            DO ji = 1, npti
348               IF    ( zGsum(ji,jl)   < rn_gstar ) THEN
349                  apartf(ji,jl) = z1_gstar * ( zGsum(ji,jl) - zGsum(ji,jl-1) ) * &
350                     &                       ( 2._wp - ( zGsum(ji,jl-1) + zGsum(ji,jl) ) * z1_gstar )
351               ELSEIF( zGsum(ji,jl-1) < rn_gstar ) THEN
352                  apartf(ji,jl) = z1_gstar * ( rn_gstar     - zGsum(ji,jl-1) ) *  &
353                     &                       ( 2._wp - ( zGsum(ji,jl-1) + rn_gstar        ) * z1_gstar )
354               ELSE
355                  apartf(ji,jl) = 0._wp
356               ENDIF
357            END DO
358         END DO
359         !
360      ELSEIF( ln_partf_exp ) THEN      !--- Exponential, more stable formulation (Lipscomb et al, 2007)
361         !                       
362         zfac = 1._wp / ( 1._wp - EXP(-z1_astar) )
363         DO jl = -1, jpl
364            DO ji = 1, npti
365               zGsum(ji,jl) = EXP( -zGsum(ji,jl) * z1_astar ) * zfac
366            END DO
367         END DO
368         DO jl = 0, jpl
369            DO ji = 1, npti
370               apartf(ji,jl) = zGsum(ji,jl-1) - zGsum(ji,jl)
371            END DO
372         END DO
373         !
374      ENDIF
375
376      !                                !--- Ridging and rafting participation concentrations
377      IF( ln_rafting .AND. ln_ridging ) THEN             !- ridging & rafting
378         DO jl = 1, jpl
379            DO ji = 1, npti
380               aridge(ji,jl) = ( 1._wp + TANH ( rn_craft * ( zhi(ji,jl) - rn_hraft ) ) ) * 0.5_wp * apartf(ji,jl)
381               araft (ji,jl) = apartf(ji,jl) - aridge(ji,jl)
382            END DO
383         END DO
384      ELSEIF( ln_ridging .AND. .NOT. ln_rafting ) THEN   !- ridging alone
385         DO jl = 1, jpl
386            DO ji = 1, npti
387               aridge(ji,jl) = apartf(ji,jl)
388               araft (ji,jl) = 0._wp
389            END DO
390         END DO
391      ELSEIF( ln_rafting .AND. .NOT. ln_ridging ) THEN   !- rafting alone   
392         DO jl = 1, jpl
393            DO ji = 1, npti
394               aridge(ji,jl) = 0._wp
395               araft (ji,jl) = apartf(ji,jl)
396            END DO
397         END DO
398      ELSE                                               !- no ridging & no rafting
399         DO jl = 1, jpl
400            DO ji = 1, npti
401               aridge(ji,jl) = 0._wp
402               araft (ji,jl) = 0._wp         
403            END DO
404         END DO
405      ENDIF
406
407      ! 2) Transfer function
408      !-----------------------------------------------------------------
409      ! Compute max and min ridged ice thickness for each ridging category.
410      ! Assume ridged ice is uniformly distributed between hrmin and hrmax.
411      !
412      ! This parameterization is a modified version of Hibler (1980).
413      ! The mean ridging thickness, zhmean, is proportional to hi^(0.5)
414      !  and for very thick ridging ice must be >= hrdg_hi_min*hi
415      !
416      ! The minimum ridging thickness, hrmin, is equal to 2*hi
417      !  (i.e., rafting) and for very thick ridging ice is
418      !  constrained by hrmin <= (zhmean + hi)/2.
419      !
420      ! The maximum ridging thickness, hrmax, is determined by zhmean and hrmin.
421      !
422      ! These modifications have the effect of reducing the ice strength
423      ! (relative to the Hibler formulation) when very thick ice is ridging.
424      !
425      ! zaksum = net area removed/ total area removed
426      ! where total area removed = area of ice that ridges
427      !         net area removed = total area removed - area of new ridges
428      !-----------------------------------------------------------------
429      zfac = 1._wp / hi_hrft
430      zaksum(1:npti) = apartf(1:npti,0)
431      !
432      DO jl = 1, jpl
433         DO ji = 1, npti
434            IF ( apartf(ji,jl) > 0._wp ) THEN
435               zhmean         = MAX( SQRT( rn_hstar * zhi(ji,jl) ), zhi(ji,jl) * hrdg_hi_min )
436               hrmin  (ji,jl) = MIN( 2._wp * zhi(ji,jl), 0.5_wp * ( zhmean + zhi(ji,jl) ) )
437               hrmax  (ji,jl) = 2._wp * zhmean - hrmin(ji,jl)
438               hraft  (ji,jl) = zhi(ji,jl) * zfac
439               hi_hrdg(ji,jl) = zhi(ji,jl) / MAX( zhmean, epsi20 )
440               !
441               ! Normalization factor : zaksum, ensures mass conservation
442               zaksum(ji) = zaksum(ji) + aridge(ji,jl) * ( 1._wp - hi_hrdg(ji,jl) )    &
443                  &                    + araft (ji,jl) * ( 1._wp - hi_hrft )
444            ELSE
445               hrmin  (ji,jl) = 0._wp 
446               hrmax  (ji,jl) = 0._wp 
447               hraft  (ji,jl) = 0._wp 
448               hi_hrdg(ji,jl) = 1._wp
449            ENDIF
450         END DO
451      END DO
452      !
453      ! 3) closing_gross
454      !-----------------
455      ! Based on the ITD of ridging and ridged ice, convert the net closing rate to a gross closing rate. 
456      ! NOTE: 0 < aksum <= 1
457      WHERE( zaksum(1:npti) > 0._wp )   ;   closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti)
458      ELSEWHERE                         ;   closing_gross(1:npti) = 0._wp
459      END WHERE
460     
461      ! correction to closing rate if excessive ice removal
462      !----------------------------------------------------
463      ! Reduce the closing rate if more than 100% of any ice category would be removed
464      ! Reduce the opening rate in proportion
465      DO jl = 1, jpl
466         DO ji = 1, npti
467            zfac = apartf(ji,jl) * closing_gross(ji) * rdt_ice
468            IF( zfac > pa_i(ji,jl) ) THEN
469               closing_gross(ji) = pa_i(ji,jl) / apartf(ji,jl) * r1_rdtice
470            ENDIF
471         END DO
472      END DO     
473
474      ! 4) correction to opening if excessive open water removal
475      !---------------------------------------------------------
476      ! Reduce the closing rate if more than 100% of the open water would be removed
477      ! Reduce the opening rate in proportion
478      DO ji = 1, npti 
479         zfac = pato_i(ji) + ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * rdt_ice
480         IF( zfac < 0._wp ) THEN           ! would lead to negative ato_i
481            opning(ji) = apartf(ji,0) * closing_gross(ji) - pato_i(ji) * r1_rdtice 
482         ELSEIF( zfac > zasum(ji) ) THEN   ! would lead to ato_i > asum
483            opning(ji) = apartf(ji,0) * closing_gross(ji) + ( zasum(ji) - pato_i(ji) ) * r1_rdtice 
484         ENDIF
485      END DO
486      !
487   END SUBROUTINE rdgrft_prep
488
489
490   SUBROUTINE rdgrft_shift
491      !!-------------------------------------------------------------------
492      !!                ***  ROUTINE rdgrft_shift ***
493      !!
494      !! ** Purpose :   shift ridging ice among thickness categories of ice thickness
495      !!
496      !! ** Method  :   Remove area, volume, and energy from each ridging category
497      !!                and add to thicker ice categories.
498      !!-------------------------------------------------------------------
499      !
500      INTEGER  ::   ji, jj, jl, jl1, jl2, jk   ! dummy loop indices
501      REAL(wp) ::   hL, hR, farea              ! left and right limits of integration and new area going to jl2
502      REAL(wp) ::   vsw                        ! vol of water trapped into ridges
503      REAL(wp) ::   afrdg, afrft               ! fraction of category area ridged/rafted
504      REAL(wp)                  ::   airdg1, oirdg1, aprdg1, virdg1, sirdg1
505      REAL(wp)                  ::   airft1, oirft1, aprft1
506      REAL(wp), DIMENSION(jpij) ::   airdg2, oirdg2, aprdg2, virdg2, sirdg2, vsrdg, vprdg  ! area etc of new ridges
507      REAL(wp), DIMENSION(jpij) ::   airft2, oirft2, aprft2, virft , sirft , vsrft, vprft  ! area etc of rafted ice
508      !
509      REAL(wp), DIMENSION(jpij) ::   ersw             ! enth of water trapped into ridges
510      REAL(wp), DIMENSION(jpij) ::   zswitch, fvol    ! new ridge volume going to jl2
511      REAL(wp), DIMENSION(jpij) ::   z1_ai            ! 1 / a
512      !
513      REAL(wp), DIMENSION(jpij,nlay_s) ::   esrft     ! snow energy of rafting ice
514      REAL(wp), DIMENSION(jpij,nlay_i) ::   eirft     ! ice  energy of rafting ice
515      REAL(wp), DIMENSION(jpij,nlay_s) ::   esrdg     ! enth*volume of new ridges     
516      REAL(wp), DIMENSION(jpij,nlay_i) ::   eirdg     ! enth*volume of new ridges
517      !
518      INTEGER , DIMENSION(jpij) ::   itest_rdg, itest_rft   ! test for conservation
519      !!-------------------------------------------------------------------
520
521      ! 1) Change in open water area due to closing and opening
522      !--------------------------------------------------------
523      DO ji = 1, npti
524         ato_i_1d(ji) = MAX( 0._wp, ato_i_1d(ji) + ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * rdt_ice )
525      END DO
526     
527      ! 2) compute categories in which ice is removed (jl1)
528      !----------------------------------------------------
529      DO jl1 = 1, jpl
530
531         CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d(1:npti), s_i(:,:,jl1) )
532
533         DO ji = 1, npti
534
535            IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN   ! only if ice is ridging
536
537               z1_ai(ji) = 1._wp / a_i_2d(ji,jl1)
538
539               ! area of ridging / rafting ice (airdg1) and of new ridge (airdg2)
540               airdg1 = aridge(ji,jl1) * closing_gross(ji) * rdt_ice
541               airft1 = araft (ji,jl1) * closing_gross(ji) * rdt_ice
542
543               airdg2(ji) = airdg1 * hi_hrdg(ji,jl1)
544               airft2(ji) = airft1 * hi_hrft
545
546               ! ridging /rafting fractions
547               afrdg = airdg1 * z1_ai(ji)
548               afrft = airft1 * z1_ai(ji)
549
550               ! volume and enthalpy (J/m2, >0) of seawater trapped into ridges
551               vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg
552               ersw(ji) = -rhoi * vsw * rcp * sst_1d(ji)   ! clem: if sst>0, then ersw <0 (is that possible?)
553
554               ! volume etc of ridging / rafting ice and new ridges (vi, vs, sm, oi, es, ei)
555               virdg1     = v_i_2d (ji,jl1)   * afrdg
556               virdg2(ji) = v_i_2d (ji,jl1)   * afrdg * ( 1. + rn_porordg )
557               vsrdg(ji)  = v_s_2d (ji,jl1)   * afrdg
558               sirdg1     = sv_i_2d(ji,jl1)   * afrdg
559               sirdg2(ji) = sv_i_2d(ji,jl1)   * afrdg + vsw * sss_1d(ji)
560               oirdg1     = oa_i_2d(ji,jl1)   * afrdg
561               oirdg2(ji) = oa_i_2d(ji,jl1)   * afrdg * hi_hrdg(ji,jl1) 
562
563               virft(ji)  = v_i_2d (ji,jl1)   * afrft
564               vsrft(ji)  = v_s_2d (ji,jl1)   * afrft
565               sirft(ji)  = sv_i_2d(ji,jl1)   * afrft 
566               oirft1     = oa_i_2d(ji,jl1)   * afrft 
567               oirft2(ji) = oa_i_2d(ji,jl1)   * afrft * hi_hrft 
568
569               IF ( ln_pnd_H12 ) THEN
570                  aprdg1     = a_ip_2d(ji,jl1) * afrdg
571                  aprdg2(ji) = a_ip_2d(ji,jl1) * afrdg * hi_hrdg(ji,jl1)
572                  vprdg (ji) = v_ip_2d(ji,jl1) * afrdg
573                  aprft1     = a_ip_2d(ji,jl1) * afrft
574                  aprft2(ji) = a_ip_2d(ji,jl1) * afrft * hi_hrft
575                  vprft (ji) = v_ip_2d(ji,jl1) * afrft
576               ENDIF
577
578               ! Ice-ocean exchanges associated with ice porosity
579               wfx_dyn_1d(ji) = wfx_dyn_1d(ji) - vsw * rhoi * r1_rdtice   ! increase in ice volume due to seawater frozen in voids
580               sfx_dyn_1d(ji) = sfx_dyn_1d(ji) - vsw * sss_1d(ji) * rhoi * r1_rdtice
581               hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ersw(ji) * r1_rdtice          ! > 0 [W.m-2]
582
583               ! Put the snow lost by ridging into the ocean
584               !  Note that esrdg > 0; the ocean must cool to melt snow. If the ocean temp = Tf already, new ice must grow.
585               wfx_snw_dyn_1d(ji) = wfx_snw_dyn_1d(ji) + ( rhos * vsrdg(ji) * ( 1._wp - rn_fsnwrdg )   &   ! fresh water source for ocean
586                  &                                      + rhos * vsrft(ji) * ( 1._wp - rn_fsnwrft ) ) * r1_rdtice
587
588               ! virtual salt flux to keep salinity constant
589               IF( nn_icesal /= 2 )  THEN
590                  sirdg2(ji)     = sirdg2(ji)     - vsw * ( sss_1d(ji) - s_i_1d(ji) )        ! ridge salinity = s_i
591                  sfx_bri_1d(ji) = sfx_bri_1d(ji) + sss_1d(ji) * vsw * rhoi * r1_rdtice  &  ! put back sss_m into the ocean
592                     &                            - s_i_1d(ji) * vsw * rhoi * r1_rdtice     ! and get  s_i  from the ocean
593               ENDIF
594
595               ! Remove area, volume of new ridge to each category jl1
596               !------------------------------------------------------
597               a_i_2d (ji,jl1) = a_i_2d (ji,jl1) - airdg1    - airft1
598               v_i_2d (ji,jl1) = v_i_2d (ji,jl1) - virdg1    - virft(ji)
599               v_s_2d (ji,jl1) = v_s_2d (ji,jl1) - vsrdg(ji) - vsrft(ji)
600               sv_i_2d(ji,jl1) = sv_i_2d(ji,jl1) - sirdg1    - sirft(ji)
601               oa_i_2d(ji,jl1) = oa_i_2d(ji,jl1) - oirdg1    - oirft1
602               IF ( ln_pnd_H12 ) THEN
603                  a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - aprdg1    - aprft1
604                  v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - vprdg(ji) - vprft(ji)
605               ENDIF
606            ENDIF
607
608         END DO ! ji
609
610         ! special loop for e_s because of layers jk
611         DO jk = 1, nlay_s
612            DO ji = 1, npti
613               IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
614                  ! Compute ridging /rafting fractions
615                  afrdg = aridge(ji,jl1) * closing_gross(ji) * rdt_ice * z1_ai(ji)
616                  afrft = araft (ji,jl1) * closing_gross(ji) * rdt_ice * z1_ai(ji)
617                  ! Compute ridging /rafting ice and new ridges for es
618                  esrdg(ji,jk) = ze_s_2d (ji,jk,jl1) * afrdg
619                  esrft(ji,jk) = ze_s_2d (ji,jk,jl1) * afrft
620                  ! Put the snow lost by ridging into the ocean
621                  hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ( - esrdg(ji,jk) * ( 1._wp - rn_fsnwrdg )   &                 ! heat sink for ocean (<0, W.m-2)
622                     &                                - esrft(ji,jk) * ( 1._wp - rn_fsnwrft ) ) * r1_rdtice
623                  !
624                  ! Remove energy of new ridge to each category jl1
625                  !-------------------------------------------------
626                  ze_s_2d(ji,jk,jl1) = ze_s_2d(ji,jk,jl1) * ( 1._wp - afrdg - afrft ) 
627               ENDIF
628            END DO
629         END DO
630                 
631         ! special loop for e_i because of layers jk
632         DO jk = 1, nlay_i
633            DO ji = 1, npti
634               IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
635                  ! Compute ridging /rafting fractions
636                  afrdg = aridge(ji,jl1) * closing_gross(ji) * rdt_ice * z1_ai(ji)
637                  afrft = araft (ji,jl1) * closing_gross(ji) * rdt_ice * z1_ai(ji)
638                  ! Compute ridging ice and new ridges for ei
639                  eirdg(ji,jk) = ze_i_2d (ji,jk,jl1) * afrdg + ersw(ji) * r1_nlay_i
640                  eirft(ji,jk) = ze_i_2d (ji,jk,jl1) * afrft
641                  !
642                  ! Remove energy of new ridge to each category jl1
643                  !-------------------------------------------------
644                  ze_i_2d(ji,jk,jl1) = ze_i_2d(ji,jk,jl1) * ( 1._wp - afrdg - afrft ) 
645               ENDIF
646            END DO
647         END DO
648         
649         ! 3) compute categories in which ice is added (jl2)
650         !--------------------------------------------------
651         itest_rdg(1:npti) = 0
652         itest_rft(1:npti) = 0
653         DO jl2  = 1, jpl 
654            !
655            DO ji = 1, npti
656
657               IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
658
659                  ! Compute the fraction of ridged ice area and volume going to thickness category jl2
660                  IF( hrmin(ji,jl1) <= hi_max(jl2) .AND. hrmax(ji,jl1) > hi_max(jl2-1) ) THEN
661                     hL = MAX( hrmin(ji,jl1), hi_max(jl2-1) )
662                     hR = MIN( hrmax(ji,jl1), hi_max(jl2)   )
663                     farea    = ( hR      - hL      ) / ( hrmax(ji,jl1)                 - hrmin(ji,jl1)                 )
664                     fvol(ji) = ( hR * hR - hL * hL ) / ( hrmax(ji,jl1) * hrmax(ji,jl1) - hrmin(ji,jl1) * hrmin(ji,jl1) )
665                     !
666                     itest_rdg(ji) = 1   ! test for conservation
667                  ELSE
668                     farea    = 0._wp 
669                     fvol(ji) = 0._wp                 
670                  ENDIF
671
672                  ! Compute the fraction of rafted ice area and volume going to thickness category jl2
673                  IF( hraft(ji,jl1) <= hi_max(jl2) .AND. hraft(ji,jl1) >  hi_max(jl2-1) ) THEN
674                     zswitch(ji) = 1._wp
675                     !
676                     itest_rft(ji) = 1   ! test for conservation
677                  ELSE
678                     zswitch(ji) = 0._wp
679                  ENDIF
680                  !
681                  ! Patch to ensure perfect conservation if ice thickness goes mad
682                  ! Sometimes thickness is larger than hi_max(jpl) because of advection scheme (for very small areas)
683                  ! Then ice volume is removed from one category but the ridging/rafting scheme
684                  ! does not know where to move it, leading to a conservation issue. 
685                  IF( itest_rdg(ji) == 0 .AND. jl2 == jpl ) THEN   ;   farea = 1._wp   ;   fvol(ji) = 1._wp   ;   ENDIF
686                  IF( itest_rft(ji) == 0 .AND. jl2 == jpl )      zswitch(ji) = 1._wp
687                  !
688                  ! Add area, volume of new ridge to category jl2
689                  !----------------------------------------------
690                  a_i_2d (ji,jl2) = a_i_2d (ji,jl2) + ( airdg2(ji) * farea    + airft2(ji) * zswitch(ji) )
691                  oa_i_2d(ji,jl2) = oa_i_2d(ji,jl2) + ( oirdg2(ji) * farea    + oirft2(ji) * zswitch(ji) )
692                  v_i_2d (ji,jl2) = v_i_2d (ji,jl2) + ( virdg2(ji) * fvol(ji) + virft (ji) * zswitch(ji) )
693                  sv_i_2d(ji,jl2) = sv_i_2d(ji,jl2) + ( sirdg2(ji) * fvol(ji) + sirft (ji) * zswitch(ji) )
694                  v_s_2d (ji,jl2) = v_s_2d (ji,jl2) + ( vsrdg (ji) * rn_fsnwrdg * fvol(ji)  +  &
695                     &                                  vsrft (ji) * rn_fsnwrft * zswitch(ji) )
696                  IF ( ln_pnd_H12 ) THEN
697                     v_ip_2d (ji,jl2) = v_ip_2d(ji,jl2) + (   vprdg (ji) * rn_fpndrdg * fvol   (ji)   &
698                        &                                   + vprft (ji) * rn_fpndrft * zswitch(ji)   )
699                     a_ip_2d (ji,jl2) = a_ip_2d(ji,jl2) + (   aprdg2(ji) * rn_fpndrdg * farea         & 
700                        &                                   + aprft2(ji) * rn_fpndrft * zswitch(ji)   )
701                  ENDIF
702                 
703               ENDIF
704
705            END DO
706            ! Add snow energy of new ridge to category jl2
707            !---------------------------------------------
708            DO jk = 1, nlay_s
709               DO ji = 1, npti
710                  IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp )   &
711                     &   ze_s_2d(ji,jk,jl2) = ze_s_2d(ji,jk,jl2) + ( esrdg(ji,jk) * rn_fsnwrdg * fvol(ji)  +  &
712                     &                                               esrft(ji,jk) * rn_fsnwrft * zswitch(ji) )
713               END DO
714            END DO
715            ! Add ice energy of new ridge to category jl2
716            !--------------------------------------------
717            DO jk = 1, nlay_i
718               DO ji = 1, npti
719                  IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp )   &
720                     &   ze_i_2d(ji,jk,jl2) = ze_i_2d(ji,jk,jl2) + eirdg(ji,jk) * fvol(ji) + eirft(ji,jk) * zswitch(ji)                 
721               END DO
722            END DO
723            !
724         END DO ! jl2
725         !
726      END DO ! jl1
727      !
728      ! In case ridging/rafting lead to very small negative values (sometimes it happens)
729      WHERE( a_i_2d(1:npti,:) < 0._wp )   a_i_2d(1:npti,:) = 0._wp
730      WHERE( v_i_2d(1:npti,:) < 0._wp )   v_i_2d(1:npti,:) = 0._wp
731      !
732   END SUBROUTINE rdgrft_shift
733
734
735   SUBROUTINE ice_strength
736      !!----------------------------------------------------------------------
737      !!                ***  ROUTINE ice_strength ***
738      !!
739      !! ** Purpose :   computes ice strength used in dynamics routines of ice thickness
740      !!
741      !! ** Method  :   Compute the strength of the ice pack, defined as the energy (J m-2)
742      !!              dissipated per unit area removed from the ice pack under compression,
743      !!              and assumed proportional to the change in potential energy caused
744      !!              by ridging. Note that only Hibler's formulation is stable and that
745      !!              ice strength has to be smoothed
746      !!----------------------------------------------------------------------
747      INTEGER             ::   ji, jj, jl  ! dummy loop indices
748      INTEGER             ::   ismooth     ! smoothing the resistance to deformation
749      INTEGER             ::   itframe     ! number of time steps for the P smoothing
750      REAL(wp)            ::   zp, z1_3    ! local scalars
751      REAL(wp), DIMENSION(jpi,jpj) ::   zworka           ! temporary array used here
752      REAL(wp), DIMENSION(jpi,jpj) ::   zstrp1, zstrp2   ! strength at previous time steps
753      !!----------------------------------------------------------------------
754      !                              !--------------------------------------------------!
755      IF( ln_str_H79 ) THEN          ! Ice strength => Hibler (1979) method             !
756      !                              !--------------------------------------------------!
757         strength(:,:) = rn_pstar * SUM( v_i(:,:,:), dim=3 ) * EXP( -rn_crhg * ( 1._wp - SUM( a_i(:,:,:), dim=3 ) ) )
758         ismooth = 1
759         !                           !--------------------------------------------------!
760      ELSE                           ! Zero strength                                    !
761         !                           !--------------------------------------------------!
762         strength(:,:) = 0._wp
763         ismooth = 0
764      ENDIF
765      !                              !--------------------------------------------------!
766      SELECT CASE( ismooth )         ! Smoothing ice strength                           !
767      !                              !--------------------------------------------------!
768      CASE( 1 )               !--- Spatial smoothing
769         DO jj = 2, jpjm1
770            DO ji = 2, jpim1
771               IF ( SUM( a_i(ji,jj,:) ) > 0._wp ) THEN
772                  zworka(ji,jj) = ( 4.0 * strength(ji,jj)              &
773                     &                  + strength(ji-1,jj) * tmask(ji-1,jj,1) + strength(ji+1,jj) * tmask(ji+1,jj,1) & 
774                     &                  + strength(ji,jj-1) * tmask(ji,jj-1,1) + strength(ji,jj+1) * tmask(ji,jj+1,1) &
775                     &            ) / ( 4.0 + tmask(ji-1,jj,1) + tmask(ji+1,jj,1) + tmask(ji,jj-1,1) + tmask(ji,jj+1,1) )
776               ELSE
777                  zworka(ji,jj) = 0._wp
778               ENDIF
779            END DO
780         END DO
781         
782         DO jj = 2, jpjm1
783            DO ji = 2, jpim1
784               strength(ji,jj) = zworka(ji,jj)
785            END DO
786         END DO
787         CALL lbc_lnk( 'icedyn_rdgrft', strength, 'T', 1. )
788         !
789      CASE( 2 )               !--- Temporal smoothing
790         IF ( kt_ice == nit000 ) THEN
791            zstrp1(:,:) = 0._wp
792            zstrp2(:,:) = 0._wp
793         ENDIF
794         !
795         DO jj = 2, jpjm1
796            DO ji = 2, jpim1
797               IF ( SUM( a_i(ji,jj,:) ) > 0._wp ) THEN
798                  itframe = 1 ! number of time steps for the running mean
799                  IF ( zstrp1(ji,jj) > 0._wp ) itframe = itframe + 1
800                  IF ( zstrp2(ji,jj) > 0._wp ) itframe = itframe + 1
801                  zp = ( strength(ji,jj) + zstrp1(ji,jj) + zstrp2(ji,jj) ) / itframe
802                  zstrp2  (ji,jj) = zstrp1  (ji,jj)
803                  zstrp1  (ji,jj) = strength(ji,jj)
804                  strength(ji,jj) = zp
805               ENDIF
806            END DO
807         END DO
808         CALL lbc_lnk( 'icedyn_rdgrft', strength, 'T', 1. )
809         !
810      END SELECT
811      !
812   END SUBROUTINE ice_strength
813
814   
815   SUBROUTINE ice_dyn_1d2d( kn )
816      !!-----------------------------------------------------------------------
817      !!                   ***  ROUTINE ice_dyn_1d2d ***
818      !!                 
819      !! ** Purpose :   move arrays from 1d to 2d and the reverse
820      !!-----------------------------------------------------------------------
821      INTEGER, INTENT(in) ::   kn   ! 1= from 2D to 1D   ;   2= from 1D to 2D
822      !
823      INTEGER ::   jl, jk   ! dummy loop indices
824      !!-----------------------------------------------------------------------
825      !
826      SELECT CASE( kn )
827      !                    !---------------------!
828      CASE( 1 )            !==  from 2D to 1D  ==!
829         !                 !---------------------!
830         ! fields used but not modified
831         CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d(1:npti), sss_m(:,:) )
832         CALL tab_2d_1d( npti, nptidx(1:npti), sst_1d(1:npti), sst_m(:,:) )
833         ! the following fields are modified in this routine
834         !!CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti), ato_i(:,:) )
835         !!CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d(1:npti,1:jpl), a_i(:,:,:) )
836         !!CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d  (1:npti,1:jpl), v_i  (:,:,:) )
837         CALL tab_3d_2d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s (:,:,:) )
838         CALL tab_3d_2d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
839         CALL tab_3d_2d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i(:,:,:) )
840         CALL tab_3d_2d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) )
841         CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) )
842         DO jl = 1, jpl
843            DO jk = 1, nlay_s
844               CALL tab_2d_1d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) )
845            END DO
846            DO jk = 1, nlay_i
847               CALL tab_2d_1d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
848            END DO
849         END DO
850         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_dyn_1d    (1:npti), sfx_dyn    (:,:) )
851         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_bri_1d    (1:npti), sfx_bri    (:,:) )
852         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_dyn_1d    (1:npti), wfx_dyn    (:,:) )
853         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_dyn_1d    (1:npti), hfx_dyn    (:,:) )
854         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) )
855         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d    (1:npti), wfx_pnd    (:,:) )
856
857         !                 !---------------------!
858      CASE( 2 )            !==  from 1D to 2D  ==!
859         !                 !---------------------!
860         CALL tab_1d_2d( npti, nptidx(1:npti), ato_i_1d(1:npti), ato_i(:,:) )
861         CALL tab_2d_3d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i (:,:,:) )
862         CALL tab_2d_3d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) )
863         CALL tab_2d_3d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s (:,:,:) )
864         CALL tab_2d_3d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
865         CALL tab_2d_3d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i(:,:,:) )
866         CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) )
867         CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) )
868         DO jl = 1, jpl
869            DO jk = 1, nlay_s
870               CALL tab_1d_2d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) )
871            END DO
872            DO jk = 1, nlay_i
873               CALL tab_1d_2d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
874            END DO
875         END DO
876         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_dyn_1d    (1:npti), sfx_dyn    (:,:) )
877         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_bri_1d    (1:npti), sfx_bri    (:,:) )
878         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_dyn_1d    (1:npti), wfx_dyn    (:,:) )
879         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_dyn_1d    (1:npti), hfx_dyn    (:,:) )
880         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) )
881         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d    (1:npti), wfx_pnd    (:,:) )
882         !
883      END SELECT
884      !
885   END SUBROUTINE ice_dyn_1d2d
886   
887
888   SUBROUTINE ice_dyn_rdgrft_init
889      !!-------------------------------------------------------------------
890      !!                  ***  ROUTINE ice_dyn_rdgrft_init ***
891      !!
892      !! ** Purpose :   Physical constants and parameters linked
893      !!                to the mechanical ice redistribution
894      !!
895      !! ** Method  :   Read the namdyn_rdgrft namelist
896      !!                and check the parameters values
897      !!                called at the first timestep (nit000)
898      !!
899      !! ** input   :   Namelist namdyn_rdgrft
900      !!-------------------------------------------------------------------
901      INTEGER :: ios                 ! Local integer output status for namelist read
902      !!
903      NAMELIST/namdyn_rdgrft/ ln_str_H79, rn_pstar, rn_crhg, &
904         &                    rn_csrdg  ,                    &
905         &                    ln_partf_lin, rn_gstar,        &
906         &                    ln_partf_exp, rn_astar,        & 
907         &                    ln_ridging, rn_hstar, rn_porordg, rn_fsnwrdg, rn_fpndrdg,  & 
908         &                    ln_rafting, rn_hraft, rn_craft  , rn_fsnwrft, rn_fpndrft
909      !!-------------------------------------------------------------------
910      !
911      REWIND( numnam_ice_ref )              ! Namelist namicetdme in reference namelist : Ice mechanical ice redistribution
912      READ  ( numnam_ice_ref, namdyn_rdgrft, IOSTAT = ios, ERR = 901)
913901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdyn_rdgrft in reference namelist', lwp )
914      REWIND( numnam_ice_cfg )              ! Namelist namdyn_rdgrft in configuration namelist : Ice mechanical ice redistribution
915      READ  ( numnam_ice_cfg, namdyn_rdgrft, IOSTAT = ios, ERR = 902 )
916902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdyn_rdgrft in configuration namelist', lwp )
917      IF(lwm) WRITE ( numoni, namdyn_rdgrft )
918      !
919      IF (lwp) THEN                          ! control print
920         WRITE(numout,*)
921         WRITE(numout,*) 'ice_dyn_rdgrft_init: ice parameters for ridging/rafting '
922         WRITE(numout,*) '~~~~~~~~~~~~~~~~~~'
923         WRITE(numout,*) '   Namelist namdyn_rdgrft:'
924         WRITE(numout,*) '      ice strength parameterization Hibler (1979)              ln_str_H79   = ', ln_str_H79 
925         WRITE(numout,*) '            1st bulk-rheology parameter                        rn_pstar     = ', rn_pstar
926         WRITE(numout,*) '            2nd bulk-rhelogy parameter                         rn_crhg      = ', rn_crhg
927         WRITE(numout,*) '      Fraction of shear energy contributing to ridging         rn_csrdg     = ', rn_csrdg 
928         WRITE(numout,*) '      linear ridging participation function                    ln_partf_lin = ', ln_partf_lin
929         WRITE(numout,*) '            Fraction of ice coverage contributing to ridging   rn_gstar     = ', rn_gstar
930         WRITE(numout,*) '      Exponential ridging participation function               ln_partf_exp = ', ln_partf_exp
931         WRITE(numout,*) '            Equivalent to G* for an exponential function       rn_astar     = ', rn_astar
932         WRITE(numout,*) '      Ridging of ice sheets or not                             ln_ridging   = ', ln_ridging
933         WRITE(numout,*) '            max ridged ice thickness                           rn_hstar     = ', rn_hstar
934         WRITE(numout,*) '            Initial porosity of ridges                         rn_porordg   = ', rn_porordg
935         WRITE(numout,*) '            Fraction of snow volume conserved during ridging   rn_fsnwrdg   = ', rn_fsnwrdg 
936         WRITE(numout,*) '            Fraction of pond volume conserved during ridging   rn_fpndrdg   = ', rn_fpndrdg 
937         WRITE(numout,*) '      Rafting of ice sheets or not                             ln_rafting   = ', ln_rafting
938         WRITE(numout,*) '            Parmeter thickness (threshold between ridge-raft)  rn_hraft     = ', rn_hraft
939         WRITE(numout,*) '            Rafting hyperbolic tangent coefficient             rn_craft     = ', rn_craft 
940         WRITE(numout,*) '            Fraction of snow volume conserved during rafting   rn_fsnwrft   = ', rn_fsnwrft 
941         WRITE(numout,*) '            Fraction of pond volume conserved during rafting   rn_fpndrft   = ', rn_fpndrft 
942      ENDIF
943      !
944      IF ( ( ln_partf_lin .AND. ln_partf_exp ) .OR. ( .NOT.ln_partf_lin .AND. .NOT.ln_partf_exp ) ) THEN
945         CALL ctl_stop( 'ice_dyn_rdgrft_init: choose one and only one participation function (ln_partf_lin or ln_partf_exp)' )
946      ENDIF
947      !                              ! allocate tke arrays
948      IF( ice_dyn_rdgrft_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'ice_dyn_rdgrft_init: unable to allocate arrays' )
949      !
950  END SUBROUTINE ice_dyn_rdgrft_init
951
952#else
953   !!----------------------------------------------------------------------
954   !!   Default option         Empty module           NO SI3 sea-ice model
955   !!----------------------------------------------------------------------
956#endif
957
958   !!======================================================================
959END MODULE icedyn_rdgrft
Note: See TracBrowser for help on using the repository browser.