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/2020/temporary_r4_trunk/src/ICE – NEMO

source: NEMO/branches/2020/temporary_r4_trunk/src/ICE/icedyn_rdgrft.F90 @ 13469

Last change on this file since 13469 was 13469, checked in by smasson, 4 years ago

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