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_r11842_SI3-10_EAP/src/ICE – NEMO

source: NEMO/branches/2019/dev_r11842_SI3-10_EAP/src/ICE/icedyn_rdgrft.F90 @ 13662

Last change on this file since 13662 was 13662, checked in by clem, 4 years ago

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