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/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/ICE – NEMO

source: NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/ICE/icedyn_rdgrft.F90 @ 14029

Last change on this file since 14029 was 14029, checked in by rlod, 4 years ago

Caught up with trunk rev 14029

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