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 branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rdgrft.F90 @ 8586

Last change on this file since 8586 was 8586, checked in by gm, 6 years ago

#1911 (ENHANCE-09): PART I.3 - phasing with branch dev_r8183_ICEMODEL revision 8575

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