source: NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icedyn_rdgrft.F90 @ 12636

Last change on this file since 12636 was 12636, checked in by dancopsey, 8 months ago

Add:

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