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_merge_2017/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rdgrft.F90 @ 9271

Last change on this file since 9271 was 9271, checked in by clem, 6 years ago

first steps for having more than 1 snow layers in the ice (in theory). There is still icethd_dh.F90 routine to change

File size: 51.1 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( ln_timing    )   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         CALL tab_3d_2d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) )
256         CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) )
257         DO jl = 1, jpl
258            DO jk = 1, nlay_s
259               CALL tab_2d_1d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) )
260            END DO
261            DO jk = 1, nlay_i
262               CALL tab_2d_1d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
263            END DO
264         END DO
265         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_dyn_1d    (1:npti), sfx_dyn    (:,:) )
266         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_bri_1d    (1:npti), sfx_bri    (:,:) )
267         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_dyn_1d    (1:npti), wfx_dyn    (:,:) )
268         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_dyn_1d    (1:npti), hfx_dyn    (:,:) )
269         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) )
270         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:) )
271         
272         !-----------------------------------------------------------------------------!
273         ! 3) Ridging iteration
274         !-----------------------------------------------------------------------------!
275         iter            = 1
276         iterate_ridging = 1     
277         DO WHILE( iterate_ridging > 0 .AND. iter < jp_itermax )
278
279            CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d )
280
281            ! 3.2 Redistribute area, volume, and energy.
282            !-----------------------------------------------------------------------------!
283            CALL rdgrft_shift
284
285            ! 3.4 Do we keep on iterating?
286            !-----------------------------------------------------------------------------!
287            ! Check whether asum = 1.  If not (because the closing and opening
288            ! rates were reduced above), ridge again with new rates.
289            iterate_ridging = 0
290            DO ji = 1, npti
291               zfac = 1._wp - ( ato_i_1d(ji) + SUM( a_i_2d(ji,:) ) )
292               IF( ABS( zfac ) < epsi10 ) THEN
293                  closing_net(ji) = 0._wp
294                  opning     (ji) = 0._wp
295                  ato_i_1d   (ji) = MAX( 0._wp, 1._wp - SUM( a_i_2d(ji,:) ) )
296               ELSE
297                  iterate_ridging  = 1
298                  zdivu_adv  (ji) = zfac * r1_rdtice
299                  closing_net(ji) = MAX( 0._wp, -zdivu_adv(ji) )
300                  opning     (ji) = MAX( 0._wp,  zdivu_adv(ji) )
301               ENDIF
302            END DO
303            !
304            iter = iter + 1
305            IF( iter  >  jp_itermax )    CALL ctl_warn( 'icedyn_rdgrft: non-converging ridging scheme' )
306            !
307         END DO
308
309         !-----------
310         ! 1D <==> 2D
311         !-----------
312         CALL tab_1d_2d( npti, nptidx(1:npti), ato_i_1d(1:npti), ato_i(:,:) )
313         CALL tab_2d_3d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i (:,:,:) )
314         CALL tab_2d_3d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) )
315         CALL tab_2d_3d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s (:,:,:) )
316         CALL tab_2d_3d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
317         CALL tab_2d_3d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i(:,:,:) )
318         CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) )
319         CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) )
320         DO jl = 1, jpl
321            DO jk = 1, nlay_s
322               CALL tab_1d_2d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) )
323            END DO
324            DO jk = 1, nlay_i
325               CALL tab_1d_2d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
326            END DO
327         END DO
328         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_dyn_1d    (1:npti), sfx_dyn    (:,:) )
329         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_bri_1d    (1:npti), sfx_bri    (:,:) )
330         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_dyn_1d    (1:npti), wfx_dyn    (:,:) )
331         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_dyn_1d    (1:npti), hfx_dyn    (:,:) )
332         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) )
333         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:) )
334
335      ENDIF ! npti > 0
336   
337      CALL ice_var_agg( 1 ) 
338
339      ! controls
340      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icedyn_rdgrft', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
341      IF( ln_ctl       )   CALL ice_prt3D   ('icedyn_rdgrft')                                                             ! prints
342      IF( ln_timing    )   CALL timing_stop ('icedyn_rdgrft')                                                             ! timing
343      !
344   END SUBROUTINE ice_dyn_rdgrft
345
346
347   SUBROUTINE rdgrft_prep( pa_i, pv_i, pato_i )
348      !!-------------------------------------------------------------------
349      !!                ***  ROUTINE rdgrft_prep ***
350      !!
351      !! ** Purpose :   preparation for ridging calculations
352      !!
353      !! ** Method  :   Compute the thickness distribution of the ice and open water
354      !!              participating in ridging and of the resulting ridges.
355      !!-------------------------------------------------------------------
356      REAL(wp), DIMENSION(:)  , INTENT(in) ::   pato_i 
357      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pa_i, pv_i 
358      !!
359      INTEGER  ::   ji, jl                     ! dummy loop indices
360      REAL(wp) ::   z1_gstar, z1_astar, zhmean, zfac   ! local scalar
361      REAL(wp), DIMENSION(jpij)        ::   zasum, z1_asum, zaksum   ! sum of a_i+ato_i and reverse
362      REAL(wp), DIMENSION(jpij,jpl)    ::   zhi                      ! ice thickness
363      REAL(wp), DIMENSION(jpij,-1:jpl) ::   zGsum                    ! zGsum(n) = sum of areas in categories 0 to n
364      !--------------------------------------------------------------------
365
366      z1_gstar = 1._wp / rn_gstar
367      z1_astar = 1._wp / rn_astar
368
369      !                       ! Ice thickness needed for rafting
370      WHERE( pa_i(1:npti,:) > 0._wp )   ;   zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:)
371      ELSEWHERE                         ;   zhi(1:npti,:) = 0._wp
372      END WHERE
373
374      !-----------------------------------------------------------------
375      ! 1) Participation function: a(h) = b(h).g(h) (apartf)
376      !-----------------------------------------------------------------
377      ! Compute the participation function apartf; this is analogous to
378      ! a(h) = b(h)g(h) as defined in Thorndike et al. (1975).
379      ! area lost from category n due to ridging/closing
380      ! apartf(n)   = total area lost due to ridging/closing
381      ! assume b(h) = (2/Gstar) * (1 - G(h)/Gstar).
382      !
383      ! The expressions for apartf are found by integrating b(h)g(h) between
384      ! the category boundaries.
385      ! apartf is always >= 0 and SUM(apartf(0:jpl))=1
386      !-----------------------------------------------------------------
387      !
388      ! Compute total area of ice plus open water.
389      ! This is in general not equal to one because of divergence during transport
390      zasum(1:npti) = pato_i(1:npti) + SUM( pa_i(1:npti,:), dim=2 )
391      !
392      WHERE( zasum(1:npti) > 0._wp )   ;   z1_asum(1:npti) = 1._wp / zasum(1:npti)
393      ELSEWHERE                        ;   z1_asum(1:npti) = 0._wp
394      END WHERE
395      ! Compute cumulative thickness distribution function
396      ! Compute the cumulative thickness distribution function zGsum,
397      ! where zGsum(n) is the fractional area in categories 0 to n.
398      ! initial value (in h = 0) equals open water area
399      zGsum(1:npti,-1) = 0._wp
400      zGsum(1:npti,0 ) = pato_i(1:npti) * z1_asum(1:npti)
401      DO jl = 1, jpl
402         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)
403      END DO
404      !
405      IF( ln_partf_lin ) THEN          !--- Linear formulation (Thorndike et al., 1975)
406         DO jl = 0, jpl   
407            DO ji = 1, npti
408               IF    ( zGsum(ji,jl)   < rn_gstar ) THEN
409                  apartf(ji,jl) = z1_gstar * ( zGsum(ji,jl) - zGsum(ji,jl-1) ) * &
410                     &                       ( 2._wp - ( zGsum(ji,jl-1) + zGsum(ji,jl) ) * z1_gstar )
411               ELSEIF( zGsum(ji,jl-1) < rn_gstar ) THEN
412                  apartf(ji,jl) = z1_gstar * ( rn_gstar     - zGsum(ji,jl-1) ) *  &
413                     &                       ( 2._wp - ( zGsum(ji,jl-1) + rn_gstar        ) * z1_gstar )
414               ELSE
415                  apartf(ji,jl) = 0._wp
416               ENDIF
417            END DO
418         END DO
419         !
420      ELSEIF( ln_partf_exp ) THEN      !--- Exponential, more stable formulation (Lipscomb et al, 2007)
421         !                       
422         zfac = 1._wp / ( 1._wp - EXP(-z1_astar) )
423         DO jl = -1, jpl
424            DO ji = 1, npti
425               zGsum(ji,jl) = EXP( -zGsum(ji,jl) * z1_astar ) * zfac
426            END DO
427         END DO
428         DO jl = 0, jpl
429            DO ji = 1, npti
430               apartf(ji,jl) = zGsum(ji,jl-1) - zGsum(ji,jl)
431            END DO
432         END DO
433         !
434      ENDIF
435
436      !                                !--- Ridging and rafting participation concentrations
437      IF( ln_rafting .AND. ln_ridging ) THEN             !- ridging & rafting
438         DO jl = 1, jpl
439            DO ji = 1, npti
440               aridge(ji,jl) = ( 1._wp + TANH ( rn_craft * ( zhi(ji,jl) - rn_hraft ) ) ) * 0.5_wp * apartf(ji,jl)
441               araft (ji,jl) = apartf(ji,jl) - aridge(ji,jl)
442            END DO
443         END DO
444      ELSEIF( ln_ridging .AND. .NOT. ln_rafting ) THEN   !- ridging alone
445         DO jl = 1, jpl
446            DO ji = 1, npti
447               aridge(ji,jl) = apartf(ji,jl)
448               araft (ji,jl) = 0._wp
449            END DO
450         END DO
451      ELSEIF( ln_rafting .AND. .NOT. ln_ridging ) THEN   !- rafting alone   
452         DO jl = 1, jpl
453            DO ji = 1, npti
454               aridge(ji,jl) = 0._wp
455               araft (ji,jl) = apartf(ji,jl)
456            END DO
457         END DO
458      ELSE                                               !- no ridging & no rafting
459         DO jl = 1, jpl
460            DO ji = 1, npti
461               aridge(ji,jl) = 0._wp
462               araft (ji,jl) = 0._wp         
463            END DO
464         END DO
465      ENDIF
466
467      !-----------------------------------------------------------------
468      ! 2) Transfer function
469      !-----------------------------------------------------------------
470      ! Compute max and min ridged ice thickness for each ridging category.
471      ! Assume ridged ice is uniformly distributed between hrmin and hrmax.
472      !
473      ! This parameterization is a modified version of Hibler (1980).
474      ! The mean ridging thickness, zhmean, is proportional to hi^(0.5)
475      !  and for very thick ridging ice must be >= hrdg_hi_min*hi
476      !
477      ! The minimum ridging thickness, hrmin, is equal to 2*hi
478      !  (i.e., rafting) and for very thick ridging ice is
479      !  constrained by hrmin <= (zhmean + hi)/2.
480      !
481      ! The maximum ridging thickness, hrmax, is determined by
482      !  zhmean and hrmin.
483      !
484      ! These modifications have the effect of reducing the ice strength
485      ! (relative to the Hibler formulation) when very thick ice is
486      ! ridging.
487      !
488      ! zaksum = net area removed/ total area removed
489      ! where total area removed = area of ice that ridges
490      !         net area removed = total area removed - area of new ridges
491      !-----------------------------------------------------------------
492
493      zfac = 1._wp / hi_hrft
494      zaksum(1:npti) = apartf(1:npti,0)
495      ! Transfer function
496      DO jl = 1, jpl !all categories have a specific transfer function
497         DO ji = 1, npti
498            IF ( apartf(ji,jl) > 0._wp ) THEN
499               zhmean         = MAX( SQRT( rn_hstar * zhi(ji,jl) ), zhi(ji,jl) * hrdg_hi_min )
500               hrmin  (ji,jl) = MIN( 2._wp * zhi(ji,jl), 0.5_wp * ( zhmean + zhi(ji,jl) ) )
501               hrmax  (ji,jl) = 2._wp * zhmean - hrmin(ji,jl)
502               hraft  (ji,jl) = zhi(ji,jl) * zfac
503               hi_hrdg(ji,jl) = zhi(ji,jl) / MAX( zhmean, epsi20 )
504               !
505               ! Normalization factor : zaksum, ensures mass conservation
506               zaksum(ji) = zaksum(ji) + aridge(ji,jl) * ( 1._wp - hi_hrdg(ji,jl) )    &
507                  &                    + araft (ji,jl) * ( 1._wp - hi_hrft )
508            ELSE
509               hrmin  (ji,jl) = 0._wp 
510               hrmax  (ji,jl) = 0._wp 
511               hraft  (ji,jl) = 0._wp 
512               hi_hrdg(ji,jl) = 1._wp
513            ENDIF
514         END DO
515      END DO
516      !
517      ! 3.1 closing_gross
518      !-----------------------------------------------------------------------------!
519      ! Based on the ITD of ridging and ridged ice, convert the net
520      !  closing rate to a gross closing rate. 
521      ! NOTE: 0 < aksum <= 1
522      WHERE( zaksum(1:npti) > 0._wp )   ;   closing_gross(1:npti) = closing_net(1:npti) / zaksum(1:npti)
523      ELSEWHERE                         ;   closing_gross(1:npti) = 0._wp
524      END WHERE
525     
526      DO ji = 1, npti 
527         ! correction to closing rate and opening if closing rate is excessive
528         !---------------------------------------------------------------------
529         ! Reduce the closing rate if more than 100% of the open water
530         ! would be removed.  Reduce the opening rate proportionately.
531         zfac = ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * rdt_ice
532         IF    ( zfac < 0._wp .AND. zfac > - pato_i(ji) ) THEN                  ! would lead to negative ato_i
533            opning(ji) = apartf(ji,0) * closing_gross(ji) - pato_i(ji) * r1_rdtice 
534         ELSEIF( zfac > 0._wp .AND. zfac > ( zasum(ji) - pato_i(ji) ) ) THEN  ! would lead to ato_i > asum
535            opning(ji) = apartf(ji,0) * closing_gross(ji) + ( zasum(ji) - pato_i(ji) ) * r1_rdtice 
536         ENDIF
537      END DO
538     
539      ! correction to closing rate / opening if excessive ice removal
540      !---------------------------------------------------------------
541      ! Reduce the closing rate if more than 100% of any ice category
542      ! would be removed.  Reduce the opening rate proportionately.
543      DO jl = 1, jpl
544         DO ji = 1, npti
545            zfac = apartf(ji,jl) * closing_gross(ji) * rdt_ice
546!!            IF( zfac > pa_i(ji,jl) .AND. zfac > 0._wp ) THEN
547            IF( zfac > pa_i(ji,jl) ) THEN
548               closing_gross(ji) = closing_gross(ji) * pa_i(ji,jl) / zfac
549            ENDIF
550         END DO
551      END DO     
552      !
553   END SUBROUTINE rdgrft_prep
554
555
556   SUBROUTINE rdgrft_shift
557      !!-------------------------------------------------------------------
558      !!                ***  ROUTINE rdgrft_shift ***
559      !!
560      !! ** Purpose :   shift ridging ice among thickness categories of ice thickness
561      !!
562      !! ** Method  :   Remove area, volume, and energy from each ridging category
563      !!              and add to thicker ice categories.
564      !!-------------------------------------------------------------------
565      !
566      INTEGER  ::   ji, jj, jl, jl1, jl2, jk   ! dummy loop indices
567      REAL(wp) ::   hL, hR, farea              ! left and right limits of integration and new area going to jl2
568      REAL(wp) ::   vsw                        ! vol of water trapped into ridges
569      REAL(wp) ::   afrdg, afrft               ! fraction of category area ridged/rafted
570      REAL(wp)                  ::   airdg1, oirdg1, aprdg1, virdg1, sirdg1
571      REAL(wp)                  ::   airft1, oirft1, aprft1
572      REAL(wp), DIMENSION(jpij) ::   airdg2, oirdg2, aprdg2, virdg2, sirdg2, vsrdg, vprdg  ! area etc of new ridges
573      REAL(wp), DIMENSION(jpij) ::   airft2, oirft2, aprft2, virft , sirft , vsrft, vprft  ! area etc of rafted ice
574      !
575      REAL(wp), DIMENSION(jpij) ::   ersw             ! enth of water trapped into ridges
576      REAL(wp), DIMENSION(jpij) ::   zswitch, fvol    ! new ridge volume going to jl2
577      REAL(wp), DIMENSION(jpij) ::   z1_ai            ! 1 / a
578      !
579      REAL(wp), DIMENSION(jpij,nlay_s) ::   esrft     ! snow energy of rafting ice
580      REAL(wp), DIMENSION(jpij,nlay_i) ::   eirft     ! ice  energy of rafting ice
581      REAL(wp), DIMENSION(jpij,nlay_s) ::   esrdg     ! enth*volume of new ridges     
582      REAL(wp), DIMENSION(jpij,nlay_i) ::   eirdg     ! enth*volume of new ridges
583      !!-------------------------------------------------------------------
584
585      !-------------------------------------------------------------------------------
586      ! 1) Compute change in open water area due to closing and opening.
587      !-------------------------------------------------------------------------------
588      DO ji = 1, npti
589         ato_i_1d(ji) = MAX( 0._wp, ato_i_1d(ji) + ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * rdt_ice )
590      END DO
591     
592      !-----------------------------------------------------------------
593      ! 2) compute categories participating to ridging/rafting (jl1)
594      !-----------------------------------------------------------------
595      DO jl1 = 1, jpl
596
597         CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d(1:npti), s_i(:,:,jl1) )
598
599         DO ji = 1, npti
600
601            !------------------------------------------------
602            ! 2.1) Identify grid cells with nonzero ridging
603            !------------------------------------------------
604            IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
605
606               z1_ai(ji) = 1._wp / a_i_2d(ji,jl1)
607
608               !--------------------------------------------------------------------
609               ! 2.2) Compute area of ridging ice (airdg1) and of new ridge (airdg2)
610               !--------------------------------------------------------------------
611               airdg1 = aridge(ji,jl1) * closing_gross(ji) * rdt_ice
612               airft1 = araft (ji,jl1) * closing_gross(ji) * rdt_ice
613
614               airdg2(ji) = airdg1 * hi_hrdg(ji,jl1)
615               airft2(ji) = airft1 * hi_hrft
616
617               !---------------------------------------------------------------
618               ! 2.3) Compute ridging /rafting fractions, make sure afrdg <=1
619               !---------------------------------------------------------------
620               afrdg = airdg1 * z1_ai(ji)
621               afrft = airft1 * z1_ai(ji)
622
623               ! volume and enthalpy (J/m2, >0) of seawater trapped into ridges
624               vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg
625               ersw(ji) = -rhoic * vsw * rcp * sst_1d(ji)   ! clem: if sst>0, then ersw <0 (is that possible?)
626
627               !---------------------------------------------------------------
628               ! 2.4) Compute ridging ice and new ridges (vi, vs, sm, oi, es, ei)
629               !---------------------------------------------------------------
630               virdg1     = v_i_2d (ji,jl1)   * afrdg
631               virdg2(ji) = v_i_2d (ji,jl1)   * afrdg * ( 1. + rn_porordg )
632               vsrdg(ji)  = v_s_2d (ji,jl1)   * afrdg
633               sirdg1     = sv_i_2d(ji,jl1)   * afrdg
634               sirdg2(ji) = sv_i_2d(ji,jl1)   * afrdg + vsw * sss_1d(ji)
635               oirdg1     = oa_i_2d(ji,jl1)   * afrdg
636               oirdg2(ji) = oa_i_2d(ji,jl1)   * afrdg * hi_hrdg(ji,jl1) 
637
638               virft(ji)  = v_i_2d (ji,jl1)   * afrft
639               vsrft(ji)  = v_s_2d (ji,jl1)   * afrft
640               sirft(ji)  = sv_i_2d(ji,jl1)   * afrft 
641               oirft1     = oa_i_2d(ji,jl1)   * afrft 
642               oirft2(ji) = oa_i_2d(ji,jl1)   * afrft * hi_hrft 
643
644               IF ( ln_pnd_H12 ) THEN
645                  aprdg1     = a_ip_2d(ji,jl1) * afrdg
646                  aprdg2(ji) = a_ip_2d(ji,jl1) * afrdg * hi_hrdg(ji,jl1)
647                  vprdg (ji) = v_ip_2d(ji,jl1) * afrdg
648                  aprft1     = a_ip_2d(ji,jl1) * afrft
649                  aprft2(ji) = a_ip_2d(ji,jl1) * afrft * hi_hrft
650                  vprft (ji) = v_ip_2d(ji,jl1) * afrft
651               ENDIF
652
653               !-----------------------------------------------------------------
654               ! 2.5) Ice-ocean exchanges
655               !-----------------------------------------------------------------
656               wfx_dyn_1d(ji) = wfx_dyn_1d(ji) - vsw * rhoic * r1_rdtice   ! increase in ice volume due to seawater frozen in voids
657               sfx_dyn_1d(ji) = sfx_dyn_1d(ji) - vsw * sss_1d(ji) * rhoic * r1_rdtice
658               hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ersw(ji) * r1_rdtice          ! > 0 [W.m-2]
659
660               ! Put the snow lost by ridging into the ocean
661               !------------------------------------------------
662               !  Note that esrdg > 0; the ocean must cool to melt snow. If the ocean temp = Tf already, new ice must grow.
663               wfx_snw_dyn_1d(ji) = wfx_snw_dyn_1d(ji) + ( rhosn * vsrdg(ji) * ( 1._wp - rn_fsnwrdg )   &   ! fresh water source for ocean
664                  &                                      + rhosn * vsrft(ji) * ( 1._wp - rn_fsnwrft ) ) * r1_rdtice
665
666               ! Put the melt pond water into the ocean
667               !------------------------------------------           
668               ! clem: I think the following lines must be commented since there
669               !       is no net mass flux between melt ponds and the ocean (see icethd_pnd.F90 for ex.)
670               !IF ( ln_pnd_fwb ) THEN
671               !   wfx_pnd_1d(ji) = wfx_pnd_1d(ji) + ( rhofw * vprdg(ji) * ( 1._wp - rn_fpndrdg )   &        ! fresh water source for ocean
672               !      &                              + rhofw * vprft(ji) * ( 1._wp - rn_fpndrft ) ) * r1_rdtice
673               !ENDIF
674
675               ! virtual salt flux to keep salinity constant
676               !------------------------------------------           
677               IF( nn_icesal /= 2 )  THEN
678                  sirdg2(ji)     = sirdg2(ji)     - vsw * ( sss_1d(ji) - s_i_1d(ji) )        ! ridge salinity = s_i
679                  sfx_bri_1d(ji) = sfx_bri_1d(ji) + sss_1d(ji) * vsw * rhoic * r1_rdtice  &  ! put back sss_m into the ocean
680                     &                            - s_i_1d(ji) * vsw * rhoic * r1_rdtice     ! and get  s_i  from the ocean
681               ENDIF
682
683
684               !------------------------------------------           
685               ! 2.6) update jl1 (removing ridged/rafted area)
686               !------------------------------------------           
687               a_i_2d (ji,jl1) = a_i_2d (ji,jl1) - airdg1    - airft1
688               v_i_2d (ji,jl1) = v_i_2d (ji,jl1) - virdg1    - virft(ji)
689               v_s_2d (ji,jl1) = v_s_2d (ji,jl1) - vsrdg(ji) - vsrft(ji)
690               sv_i_2d(ji,jl1) = sv_i_2d(ji,jl1) - sirdg1    - sirft(ji)
691               oa_i_2d(ji,jl1) = oa_i_2d(ji,jl1) - oirdg1    - oirft1
692               IF ( ln_pnd_H12 ) THEN
693                  a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - aprdg1    - aprft1
694                  v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - vprdg(ji) - vprft(ji)
695               ENDIF
696            ENDIF
697
698         END DO ! ji
699
700         ! special loop for e_s because of layers jk
701         DO jk = 1, nlay_s
702            DO ji = 1, npti
703               IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
704                  ! Compute ridging /rafting fractions
705                  afrdg = aridge(ji,jl1) * closing_gross(ji) * rdt_ice * z1_ai(ji)
706                  afrft = araft (ji,jl1) * closing_gross(ji) * rdt_ice * z1_ai(ji)
707                  ! Compute ridging ice and new ridges for es
708                  esrdg(ji,jk) = ze_s_2d (ji,jk,jl1) * afrdg
709                  esrft(ji,jk) = ze_s_2d (ji,jk,jl1) * afrft
710                  ! Put the snow lost by ridging into the ocean
711                  hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ( - esrdg(ji,jk) * ( 1._wp - rn_fsnwrdg )   &                 ! heat sink for ocean (<0, W.m-2)
712                     &                                - esrft(ji,jk) * ( 1._wp - rn_fsnwrft ) ) * r1_rdtice
713                  ! Update jl1
714                  ze_s_2d(ji,jk,jl1) = ze_s_2d(ji,jk,jl1) * ( 1._wp - afrdg - afrft ) 
715               ENDIF
716            END DO
717         END DO
718                 
719         ! special loop for e_i because of layers jk
720         DO jk = 1, nlay_i
721            DO ji = 1, npti
722               IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
723                  ! Compute ridging /rafting fractions
724                  afrdg = aridge(ji,jl1) * closing_gross(ji) * rdt_ice * z1_ai(ji)
725                  afrft = araft (ji,jl1) * closing_gross(ji) * rdt_ice * z1_ai(ji)
726                  ! Compute ridging ice and new ridges for ei
727                  eirdg(ji,jk) = ze_i_2d (ji,jk,jl1) * afrdg + ersw(ji) * r1_nlay_i
728                  eirft(ji,jk) = ze_i_2d (ji,jk,jl1) * afrft
729                  ! Update jl1
730                  ze_i_2d(ji,jk,jl1) = ze_i_2d(ji,jk,jl1) * ( 1._wp - afrdg - afrft ) 
731               ENDIF
732            END DO
733         END DO
734         
735
736         !-------------------------------------------------------------------------------
737         ! 3) Add area, volume, and energy of new ridge to each category jl2
738         !-------------------------------------------------------------------------------
739         DO jl2  = 1, jpl 
740            !
741            DO ji = 1, npti
742
743               IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
744
745                  ! Compute the fraction of ridged ice area and volume going to thickness category jl2.
746                  IF( hrmin(ji,jl1) <= hi_max(jl2) .AND. hrmax(ji,jl1) > hi_max(jl2-1) ) THEN
747                     !-----------------------------------------------------------------
748                     ! 2.8 Compute quantities used to apportion ice among categories
749                     ! in the n2 loop below
750                     !-----------------------------------------------------------------
751                     hL = MAX( hrmin(ji,jl1), hi_max(jl2-1) )
752                     hR = MIN( hrmax(ji,jl1), hi_max(jl2)   )
753                     farea    = ( hR      - hL      ) / ( hrmax(ji,jl1)                 - hrmin(ji,jl1)                 )
754                     fvol(ji) = ( hR * hR - hL * hL ) / ( hrmax(ji,jl1) * hrmax(ji,jl1) - hrmin(ji,jl1) * hrmin(ji,jl1) )
755                  ELSE
756                     farea    = 0._wp 
757                     fvol(ji) = 0._wp                 
758                  ENDIF
759
760                  ! Compute the fraction of rafted ice area and volume going to thickness category jl2
761                  !!gm see above               IF( hraft(ji) <= hi_max(jl2) .AND. hraft(ji) >  hi_max(jl2-1) ) THEN
762                  IF( hi_max(jl2-1) < hraft(ji,jl1) .AND. hraft(ji,jl1) <= hi_max(jl2)  ) THEN   ;   zswitch(ji) = 1._wp
763                  ELSE                                                                           ;   zswitch(ji) = 0._wp
764                  ENDIF
765                  !
766                  a_i_2d (ji,jl2) = a_i_2d (ji,jl2) + ( airdg2(ji) * farea    + airft2(ji) * zswitch(ji) )
767                  oa_i_2d(ji,jl2) = oa_i_2d(ji,jl2) + ( oirdg2(ji) * farea    + oirft2(ji) * zswitch(ji) )
768                  v_i_2d (ji,jl2) = v_i_2d (ji,jl2) + ( virdg2(ji) * fvol(ji) + virft (ji) * zswitch(ji) )
769                  sv_i_2d(ji,jl2) = sv_i_2d(ji,jl2) + ( sirdg2(ji) * fvol(ji) + sirft (ji) * zswitch(ji) )
770                  v_s_2d (ji,jl2) = v_s_2d (ji,jl2) + ( vsrdg (ji) * rn_fsnwrdg * fvol(ji)  +  &
771                     &                                  vsrft (ji) * rn_fsnwrft * zswitch(ji) )
772                  IF ( ln_pnd_H12 ) THEN
773                     v_ip_2d (ji,jl2) = v_ip_2d(ji,jl2) + (   vprdg (ji) * rn_fpndrdg * fvol   (ji)   &
774                        &                                   + vprft (ji) * rn_fpndrft * zswitch(ji)   )
775                     a_ip_2d (ji,jl2) = a_ip_2d(ji,jl2) + (   aprdg2(ji) * rn_fpndrdg * farea         & 
776                        &                                   + aprft2(ji) * rn_fpndrft * zswitch(ji)   )
777                  ENDIF
778                 
779               ENDIF
780
781            END DO
782            ! for snow enthalpy
783            DO jk = 1, nlay_s
784               DO ji = 1, npti
785                  IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp )   &
786                     &   ze_s_2d(ji,jk,jl2) = ze_s_2d(ji,jk,jl2) + ( esrdg(ji,jk) * rn_fsnwrdg * fvol(ji)  +  &
787                     &                                               esrft(ji,jk) * rn_fsnwrft * zswitch(ji) )
788               END DO
789            END DO
790            ! for ice enthalpy
791            DO jk = 1, nlay_i
792               DO ji = 1, npti
793                  IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp )   &
794                     &   ze_i_2d(ji,jk,jl2) = ze_i_2d(ji,jk,jl2) + eirdg(ji,jk) * fvol(ji) + eirft(ji,jk) * zswitch(ji)                 
795               END DO
796            END DO
797            !
798         END DO ! jl2
799         !
800      END DO ! jl1
801      !
802      ! In case ridging/rafting lead to very small negative values (sometimes it happens)
803      WHERE( a_i_2d(1:npti,:) < 0._wp )   a_i_2d(1:npti,:) = 0._wp
804      WHERE( v_i_2d(1:npti,:) < 0._wp )   v_i_2d(1:npti,:) = 0._wp
805      !
806   END SUBROUTINE rdgrft_shift
807
808
809   SUBROUTINE ice_strength
810      !!----------------------------------------------------------------------
811      !!                ***  ROUTINE ice_strength ***
812      !!
813      !! ** Purpose :   computes ice strength used in dynamics routines of ice thickness
814      !!
815      !! ** Method  :   Compute the strength of the ice pack, defined as the energy (J m-2)
816      !!              dissipated per unit area removed from the ice pack under compression,
817      !!              and assumed proportional to the change in potential energy caused
818      !!              by ridging. Note that only Hibler's formulation is stable and that
819      !!              ice strength has to be smoothed
820      !!----------------------------------------------------------------------
821      INTEGER             ::   ji, jj, jl  ! dummy loop indices
822      INTEGER             ::   ismooth     ! smoothing the resistance to deformation
823      INTEGER             ::   itframe     ! number of time steps for the P smoothing
824      REAL(wp)            ::   zp, z1_3    ! local scalars
825      REAL(wp), DIMENSION(jpi,jpj) ::   zworka           ! temporary array used here
826      REAL(wp), DIMENSION(jpi,jpj) ::   zstrp1, zstrp2   ! strength at previous time steps
827      !!----------------------------------------------------------------------
828      !                              !--------------------------------------------------!
829      IF( ln_str_H79 ) THEN          ! Ice strength => Hibler (1979) method             !
830      !                              !--------------------------------------------------!
831         strength(:,:) = rn_pstar * SUM( v_i(:,:,:), dim=3 ) * EXP( -rn_crhg * ( 1._wp - SUM( a_i(:,:,:), dim=3 ) ) )
832         ismooth = 1
833         !                           !--------------------------------------------------!
834      ELSE                           ! Zero strength                                    !
835         !                           !--------------------------------------------------!
836         strength(:,:) = 0._wp
837         ismooth = 0
838      ENDIF
839      !                              !--------------------------------------------------!
840      SELECT CASE( ismooth )         ! Smoothing ice strength                           !
841      !                              !--------------------------------------------------!
842      CASE( 1 )               !--- Spatial smoothing
843         DO jj = 2, jpjm1
844            DO ji = 2, jpim1
845               IF ( SUM( a_i(ji,jj,:) ) > 0._wp ) THEN
846                  zworka(ji,jj) = ( 4.0 * strength(ji,jj)              &
847                     &                  + strength(ji-1,jj) * tmask(ji-1,jj,1) + strength(ji+1,jj) * tmask(ji+1,jj,1) & 
848                     &                  + strength(ji,jj-1) * tmask(ji,jj-1,1) + strength(ji,jj+1) * tmask(ji,jj+1,1) &
849                     &            ) / ( 4.0 + tmask(ji-1,jj,1) + tmask(ji+1,jj,1) + tmask(ji,jj-1,1) + tmask(ji,jj+1,1) )
850               ELSE
851                  zworka(ji,jj) = 0._wp
852               ENDIF
853            END DO
854         END DO
855         
856         DO jj = 2, jpjm1
857            DO ji = 2, jpim1
858               strength(ji,jj) = zworka(ji,jj)
859            END DO
860         END DO
861         CALL lbc_lnk( strength, 'T', 1. )
862         !
863      CASE( 2 )               !--- Temporal smoothing
864         IF ( kt_ice == nit000 ) THEN
865            zstrp1(:,:) = 0._wp
866            zstrp2(:,:) = 0._wp
867         ENDIF
868         !
869         DO jj = 2, jpjm1
870            DO ji = 2, jpim1
871               IF ( SUM( a_i(ji,jj,:) ) > 0._wp ) THEN
872                  itframe = 1 ! number of time steps for the running mean
873                  IF ( zstrp1(ji,jj) > 0._wp ) itframe = itframe + 1
874                  IF ( zstrp2(ji,jj) > 0._wp ) itframe = itframe + 1
875                  zp = ( strength(ji,jj) + zstrp1(ji,jj) + zstrp2(ji,jj) ) / itframe
876                  zstrp2  (ji,jj) = zstrp1  (ji,jj)
877                  zstrp1  (ji,jj) = strength(ji,jj)
878                  strength(ji,jj) = zp
879               ENDIF
880            END DO
881         END DO
882         CALL lbc_lnk( strength, 'T', 1. )
883         !
884      END SELECT
885      !
886   END SUBROUTINE ice_strength
887
888
889   SUBROUTINE ice_dyn_rdgrft_init
890      !!-------------------------------------------------------------------
891      !!                  ***  ROUTINE ice_dyn_rdgrft_init ***
892      !!
893      !! ** Purpose :   Physical constants and parameters linked
894      !!                to the mechanical ice redistribution
895      !!
896      !! ** Method  :   Read the namdyn_rdgrft namelist
897      !!                and check the parameters values
898      !!                called at the first timestep (nit000)
899      !!
900      !! ** input   :   Namelist namdyn_rdgrft
901      !!-------------------------------------------------------------------
902      INTEGER :: ios                 ! Local integer output status for namelist read
903      !!
904      NAMELIST/namdyn_rdgrft/ ln_str_H79, rn_pstar, rn_crhg, &
905         &                    rn_csrdg  ,                    &
906         &                    ln_partf_lin, rn_gstar,        &
907         &                    ln_partf_exp, rn_astar,        & 
908         &                    ln_ridging, rn_hstar, rn_porordg, rn_fsnwrdg, rn_fpndrdg,  & 
909         &                    ln_rafting, rn_hraft, rn_craft  , rn_fsnwrft, rn_fpndrft
910      !!-------------------------------------------------------------------
911      !
912      REWIND( numnam_ice_ref )              ! Namelist namicetdme in reference namelist : Ice mechanical ice redistribution
913      READ  ( numnam_ice_ref, namdyn_rdgrft, IOSTAT = ios, ERR = 901)
914901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdyn_rdgrft in reference namelist', lwp )
915      REWIND( numnam_ice_cfg )              ! Namelist namdyn_rdgrft in configuration namelist : Ice mechanical ice redistribution
916      READ  ( numnam_ice_cfg, namdyn_rdgrft, IOSTAT = ios, ERR = 902 )
917902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdyn_rdgrft in configuration namelist', lwp )
918      IF(lwm) WRITE ( numoni, namdyn_rdgrft )
919      !
920      IF (lwp) THEN                          ! control print
921         WRITE(numout,*)
922         WRITE(numout,*) 'ice_dyn_rdgrft_init: ice parameters for ridging/rafting '
923         WRITE(numout,*) '~~~~~~~~~~~~~~~~~~'
924         WRITE(numout,*) '   Namelist namdyn_rdgrft:'
925         WRITE(numout,*) '      ice strength parameterization Hibler (1979)              ln_str_H79   = ', ln_str_H79 
926         WRITE(numout,*) '            1st bulk-rheology parameter                        rn_pstar     = ', rn_pstar
927         WRITE(numout,*) '            2nd bulk-rhelogy parameter                         rn_crhg      = ', rn_crhg
928         WRITE(numout,*) '      Fraction of shear energy contributing to ridging         rn_csrdg     = ', rn_csrdg 
929         WRITE(numout,*) '      linear ridging participation function                    ln_partf_lin = ', ln_partf_lin
930         WRITE(numout,*) '            Fraction of ice coverage contributing to ridging   rn_gstar     = ', rn_gstar
931         WRITE(numout,*) '      Exponential ridging participation function               ln_partf_exp = ', ln_partf_exp
932         WRITE(numout,*) '            Equivalent to G* for an exponential function       rn_astar     = ', rn_astar
933         WRITE(numout,*) '      Ridging of ice sheets or not                             ln_ridging   = ', ln_ridging
934         WRITE(numout,*) '            max ridged ice thickness                           rn_hstar     = ', rn_hstar
935         WRITE(numout,*) '            Initial porosity of ridges                         rn_porordg   = ', rn_porordg
936         WRITE(numout,*) '            Fraction of snow volume conserved during ridging   rn_fsnwrdg   = ', rn_fsnwrdg 
937         WRITE(numout,*) '            Fraction of pond volume conserved during ridging   rn_fpndrdg   = ', rn_fpndrdg 
938         WRITE(numout,*) '      Rafting of ice sheets or not                             ln_rafting   = ', ln_rafting
939         WRITE(numout,*) '            Parmeter thickness (threshold between ridge-raft)  rn_hraft     = ', rn_hraft
940         WRITE(numout,*) '            Rafting hyperbolic tangent coefficient             rn_craft     = ', rn_craft 
941         WRITE(numout,*) '            Fraction of snow volume conserved during rafting   rn_fsnwrft   = ', rn_fsnwrft 
942         WRITE(numout,*) '            Fraction of pond volume conserved during rafting   rn_fpndrft   = ', rn_fpndrft 
943      ENDIF
944      !
945      IF ( ( ln_partf_lin .AND. ln_partf_exp ) .OR. ( .NOT.ln_partf_lin .AND. .NOT.ln_partf_exp ) ) THEN
946         CALL ctl_stop( 'ice_dyn_rdgrft_init: choose one and only one participation function (ln_partf_lin or ln_partf_exp)' )
947      ENDIF
948      !                              ! allocate tke arrays
949      IF( ice_dyn_rdgrft_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'ice_dyn_rdgrft_init: unable to allocate arrays' )
950      !
951  END SUBROUTINE ice_dyn_rdgrft_init
952
953#else
954   !!----------------------------------------------------------------------
955   !!   Default option         Empty module           NO ESIM sea-ice model
956   !!----------------------------------------------------------------------
957#endif
958
959   !!======================================================================
960END MODULE icedyn_rdgrft
Note: See TracBrowser for help on using the repository browser.