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 @ 9169

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

dev_merge_2017: all SRC: finalize the removal of useless warning when reading namelist_cfg + remove all nn_closea + nn_msh replaced by a logical

File size: 50.0 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, esrdg  ! area etc of new ridges
573      REAL(wp), DIMENSION(jpij) ::   airft2, oirft2, aprft2, virft , sirft , vsrft, vprft, esrft  ! 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_i) ::   eirft     ! ice energy of rafting ice
580      REAL(wp), DIMENSION(jpij,nlay_i) ::   eirdg     ! enth*volume of new ridges
581      !!-------------------------------------------------------------------
582
583      !-------------------------------------------------------------------------------
584      ! 1) Compute change in open water area due to closing and opening.
585      !-------------------------------------------------------------------------------
586      DO ji = 1, npti
587         ato_i_1d(ji) = MAX( 0._wp, ato_i_1d(ji) + ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * rdt_ice )
588      END DO
589     
590      !-----------------------------------------------------------------
591      ! 2) compute categories participating to ridging/rafting (jl1)
592      !-----------------------------------------------------------------
593      DO jl1 = 1, jpl
594
595         CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d(1:npti), s_i(:,:,jl1) )
596
597         DO ji = 1, npti
598
599            !------------------------------------------------
600            ! 2.1) Identify grid cells with nonzero ridging
601            !------------------------------------------------
602            IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
603
604               z1_ai(ji) = 1._wp / a_i_2d(ji,jl1)
605
606               !--------------------------------------------------------------------
607               ! 2.2) Compute area of ridging ice (airdg1) and of new ridge (airdg2)
608               !--------------------------------------------------------------------
609               airdg1 = aridge(ji,jl1) * closing_gross(ji) * rdt_ice
610               airft1 = araft (ji,jl1) * closing_gross(ji) * rdt_ice
611
612               airdg2(ji) = airdg1 * hi_hrdg(ji,jl1)
613               airft2(ji) = airft1 * hi_hrft
614
615               !---------------------------------------------------------------
616               ! 2.3) Compute ridging /rafting fractions, make sure afrdg <=1
617               !---------------------------------------------------------------
618               afrdg = airdg1 * z1_ai(ji)
619               afrft = airft1 * z1_ai(ji)
620
621               ! volume and enthalpy (J/m2, >0) of seawater trapped into ridges
622               vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg
623               ersw(ji) = -rhoic * vsw * rcp * sst_1d(ji)   ! clem: if sst>0, then ersw <0 (is that possible?)
624
625               !---------------------------------------------------------------
626               ! 2.4) Compute ridging ice and new ridges (vi, vs, sm, oi, es, ei)
627               !---------------------------------------------------------------
628               virdg1     = v_i_2d (ji,jl1)   * afrdg
629               virdg2(ji) = v_i_2d (ji,jl1)   * afrdg * ( 1. + rn_porordg )
630               vsrdg(ji)  = v_s_2d (ji,jl1)   * afrdg
631               sirdg1     = sv_i_2d(ji,jl1)   * afrdg
632               sirdg2(ji) = sv_i_2d(ji,jl1)   * afrdg + vsw * sss_1d(ji)
633               oirdg1     = oa_i_2d(ji,jl1)   * afrdg
634               oirdg2(ji) = oa_i_2d(ji,jl1)   * afrdg * hi_hrdg(ji,jl1) 
635               esrdg(ji)  = ze_s_2d(ji,1,jl1) * afrdg
636
637               virft(ji)  = v_i_2d (ji,jl1)   * afrft
638               vsrft(ji)  = v_s_2d (ji,jl1)   * afrft
639               sirft(ji)  = sv_i_2d(ji,jl1)   * afrft 
640               oirft1     = oa_i_2d(ji,jl1)   * afrft 
641               oirft2(ji) = oa_i_2d(ji,jl1)   * afrft * hi_hrft 
642               esrft(ji)  = ze_s_2d(ji,1,jl1) * afrft
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               hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ( - esrdg(ji) * ( 1._wp - rn_fsnwrdg )   &                 ! heat sink for ocean (<0, W.m-2)
667                  &                                - esrft(ji) * ( 1._wp - rn_fsnwrft ) ) * r1_rdtice
668
669               ! Put the melt pond water into the ocean
670               !------------------------------------------           
671               ! clem: I think the following lines must be commented since there
672               !       is no net mass flux between melt ponds and the ocean (see icethd_pnd.F90 for ex.)
673               !IF ( ln_pnd_fwb ) THEN
674               !   wfx_pnd_1d(ji) = wfx_pnd_1d(ji) + ( rhofw * vprdg(ji) * ( 1._wp - rn_fpndrdg )   &        ! fresh water source for ocean
675               !      &                              + rhofw * vprft(ji) * ( 1._wp - rn_fpndrft ) ) * r1_rdtice
676               !ENDIF
677
678               ! virtual salt flux to keep salinity constant
679               !------------------------------------------           
680               IF( nn_icesal /= 2 )  THEN
681                  sirdg2(ji)     = sirdg2(ji)     - vsw * ( sss_1d(ji) - s_i_1d(ji) )        ! ridge salinity = s_i
682                  sfx_bri_1d(ji) = sfx_bri_1d(ji) + sss_1d(ji) * vsw * rhoic * r1_rdtice  &  ! put back sss_m into the ocean
683                     &                            - s_i_1d(ji) * vsw * rhoic * r1_rdtice     ! and get  s_i  from the ocean
684               ENDIF
685
686
687               !------------------------------------------           
688               ! 2.6) update jl1 (removing ridged/rafted area)
689               !------------------------------------------           
690               a_i_2d (ji,jl1) = a_i_2d (ji,jl1) - airdg1    - airft1
691               v_i_2d (ji,jl1) = v_i_2d (ji,jl1) - virdg1    - virft(ji)
692               v_s_2d (ji,jl1) = v_s_2d (ji,jl1) - vsrdg(ji) - vsrft(ji)
693               sv_i_2d(ji,jl1) = sv_i_2d(ji,jl1) - sirdg1    - sirft(ji)
694               oa_i_2d(ji,jl1) = oa_i_2d(ji,jl1) - oirdg1    - oirft1
695               IF ( ln_pnd_H12 ) THEN
696                  a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - aprdg1    - aprft1
697                  v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - vprdg(ji) - vprft(ji)
698               ENDIF
699               ze_s_2d(ji,1,jl1) = ze_s_2d(ji,1,jl1) - esrdg (ji) - esrft (ji)
700            ENDIF
701
702         END DO ! ji
703
704         ! special loop for e_i because of layers jk
705         DO jk = 1, nlay_i
706            DO ji = 1, npti
707               IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
708                  ! Compute ridging /rafting fractions
709                  afrdg = aridge(ji,jl1) * closing_gross(ji) * rdt_ice * z1_ai(ji)
710                  afrft = araft (ji,jl1) * closing_gross(ji) * rdt_ice * z1_ai(ji)
711                  ! Compute ridging ice and new ridges for ei
712                  eirdg(ji,jk) = ze_i_2d (ji,jk,jl1) * afrdg + ersw(ji) * r1_nlay_i
713                  eirft(ji,jk) = ze_i_2d (ji,jk,jl1) * afrft
714                  ! Update jl1
715                  ze_i_2d(ji,jk,jl1) = ze_i_2d(ji,jk,jl1) * ( 1._wp - afrdg - afrft ) 
716               ENDIF
717            END DO
718         END DO
719         
720
721         !-------------------------------------------------------------------------------
722         ! 3) Add area, volume, and energy of new ridge to each category jl2
723         !-------------------------------------------------------------------------------
724         DO jl2  = 1, jpl 
725            !
726            DO ji = 1, npti
727
728               IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
729
730                  ! Compute the fraction of ridged ice area and volume going to thickness category jl2.
731                  IF( hrmin(ji,jl1) <= hi_max(jl2) .AND. hrmax(ji,jl1) > hi_max(jl2-1) ) THEN
732                     !-----------------------------------------------------------------
733                     ! 2.8 Compute quantities used to apportion ice among categories
734                     ! in the n2 loop below
735                     !-----------------------------------------------------------------
736                     hL = MAX( hrmin(ji,jl1), hi_max(jl2-1) )
737                     hR = MIN( hrmax(ji,jl1), hi_max(jl2)   )
738                     farea    = ( hR      - hL      ) / ( hrmax(ji,jl1)                 - hrmin(ji,jl1)                 )
739                     fvol(ji) = ( hR * hR - hL * hL ) / ( hrmax(ji,jl1) * hrmax(ji,jl1) - hrmin(ji,jl1) * hrmin(ji,jl1) )
740                  ELSE
741                     farea    = 0._wp 
742                     fvol(ji) = 0._wp                 
743                  ENDIF
744
745                  ! Compute the fraction of rafted ice area and volume going to thickness category jl2
746                  !!gm see above               IF( hraft(ji) <= hi_max(jl2) .AND. hraft(ji) >  hi_max(jl2-1) ) THEN
747                  IF( hi_max(jl2-1) < hraft(ji,jl1) .AND. hraft(ji,jl1) <= hi_max(jl2)  ) THEN   ;   zswitch(ji) = 1._wp
748                  ELSE                                                                           ;   zswitch(ji) = 0._wp
749                  ENDIF
750                  !
751                  a_i_2d (ji,jl2) = a_i_2d (ji,jl2) + ( airdg2(ji) * farea    + airft2(ji) * zswitch(ji) )
752                  oa_i_2d(ji,jl2) = oa_i_2d(ji,jl2) + ( oirdg2(ji) * farea    + oirft2(ji) * zswitch(ji) )
753                  v_i_2d (ji,jl2) = v_i_2d (ji,jl2) + ( virdg2(ji) * fvol(ji) + virft (ji) * zswitch(ji) )
754                  sv_i_2d(ji,jl2) = sv_i_2d(ji,jl2) + ( sirdg2(ji) * fvol(ji) + sirft (ji) * zswitch(ji) )
755                  v_s_2d (ji,jl2) = v_s_2d (ji,jl2) + ( vsrdg (ji) * rn_fsnwrdg * fvol(ji)  +  &
756                     &                                    vsrft (ji) * rn_fsnwrft * zswitch(ji) )
757                  IF ( ln_pnd_H12 ) THEN
758                     v_ip_2d (ji,jl2) = v_ip_2d(ji,jl2) + (   vprdg (ji) * rn_fpndrdg * fvol   (ji)   &
759                        &                                   + vprft (ji) * rn_fpndrft * zswitch(ji)   )
760                     a_ip_2d (ji,jl2) = a_ip_2d(ji,jl2) + (   aprdg2(ji) * rn_fpndrdg * farea         & 
761                        &                                   + aprft2(ji) * rn_fpndrft * zswitch(ji)   )
762                  ENDIF
763                  ze_s_2d(ji,1,jl2) = ze_s_2d(ji,1,jl2) + ( esrdg (ji) * rn_fsnwrdg * fvol(ji)  +  &
764                     &                                      esrft (ji) * rn_fsnwrft * zswitch(ji) )
765                 
766               ENDIF
767
768            END DO
769
770            DO jk = 1, nlay_i
771               DO ji = 1, npti
772                  IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp )   &
773                     &   ze_i_2d(ji,jk,jl2) = ze_i_2d(ji,jk,jl2) + eirdg(ji,jk) * fvol(ji) + eirft(ji,jk) * zswitch(ji)                 
774               END DO
775            END DO
776            !
777         END DO ! jl2
778         !
779      END DO ! jl1
780      !
781      ! In case ridging/rafting lead to very small negative values (sometimes it happens)
782      WHERE( a_i_2d(1:npti,:) < 0._wp )   a_i_2d(1:npti,:) = 0._wp
783      WHERE( v_i_2d(1:npti,:) < 0._wp )   v_i_2d(1:npti,:) = 0._wp
784      !
785   END SUBROUTINE rdgrft_shift
786
787
788   SUBROUTINE ice_strength
789      !!----------------------------------------------------------------------
790      !!                ***  ROUTINE ice_strength ***
791      !!
792      !! ** Purpose :   computes ice strength used in dynamics routines of ice thickness
793      !!
794      !! ** Method  :   Compute the strength of the ice pack, defined as the energy (J m-2)
795      !!              dissipated per unit area removed from the ice pack under compression,
796      !!              and assumed proportional to the change in potential energy caused
797      !!              by ridging. Note that only Hibler's formulation is stable and that
798      !!              ice strength has to be smoothed
799      !!----------------------------------------------------------------------
800      INTEGER             ::   ji, jj, jl  ! dummy loop indices
801      INTEGER             ::   ismooth     ! smoothing the resistance to deformation
802      INTEGER             ::   itframe     ! number of time steps for the P smoothing
803      REAL(wp)            ::   zp, z1_3    ! local scalars
804      REAL(wp), DIMENSION(jpi,jpj) ::   zworka           ! temporary array used here
805      REAL(wp), DIMENSION(jpi,jpj) ::   zstrp1, zstrp2   ! strength at previous time steps
806      !!----------------------------------------------------------------------
807      !                              !--------------------------------------------------!
808      IF( ln_str_H79 ) THEN          ! Ice strength => Hibler (1979) method             !
809      !                              !--------------------------------------------------!
810         strength(:,:) = rn_pstar * SUM( v_i(:,:,:), dim=3 ) * EXP( -rn_crhg * ( 1._wp - SUM( a_i(:,:,:), dim=3 ) ) )
811         ismooth = 1
812         !                           !--------------------------------------------------!
813      ELSE                           ! Zero strength                                    !
814         !                           !--------------------------------------------------!
815         strength(:,:) = 0._wp
816         ismooth = 0
817      ENDIF
818      !                              !--------------------------------------------------!
819      SELECT CASE( ismooth )         ! Smoothing ice strength                           !
820      !                              !--------------------------------------------------!
821      CASE( 1 )               !--- Spatial smoothing
822         DO jj = 2, jpjm1
823            DO ji = 2, jpim1
824               IF ( SUM( a_i(ji,jj,:) ) > 0._wp ) THEN
825                  zworka(ji,jj) = ( 4.0 * strength(ji,jj)              &
826                     &                  + strength(ji-1,jj) * tmask(ji-1,jj,1) + strength(ji+1,jj) * tmask(ji+1,jj,1) & 
827                     &                  + strength(ji,jj-1) * tmask(ji,jj-1,1) + strength(ji,jj+1) * tmask(ji,jj+1,1) &
828                     &            ) / ( 4.0 + tmask(ji-1,jj,1) + tmask(ji+1,jj,1) + tmask(ji,jj-1,1) + tmask(ji,jj+1,1) )
829               ELSE
830                  zworka(ji,jj) = 0._wp
831               ENDIF
832            END DO
833         END DO
834         
835         DO jj = 2, jpjm1
836            DO ji = 2, jpim1
837               strength(ji,jj) = zworka(ji,jj)
838            END DO
839         END DO
840         CALL lbc_lnk( strength, 'T', 1. )
841         !
842      CASE( 2 )               !--- Temporal smoothing
843         IF ( kt_ice == nit000 ) THEN
844            zstrp1(:,:) = 0._wp
845            zstrp2(:,:) = 0._wp
846         ENDIF
847         !
848         DO jj = 2, jpjm1
849            DO ji = 2, jpim1
850               IF ( SUM( a_i(ji,jj,:) ) > 0._wp ) THEN
851                  itframe = 1 ! number of time steps for the running mean
852                  IF ( zstrp1(ji,jj) > 0._wp ) itframe = itframe + 1
853                  IF ( zstrp2(ji,jj) > 0._wp ) itframe = itframe + 1
854                  zp = ( strength(ji,jj) + zstrp1(ji,jj) + zstrp2(ji,jj) ) / itframe
855                  zstrp2  (ji,jj) = zstrp1  (ji,jj)
856                  zstrp1  (ji,jj) = strength(ji,jj)
857                  strength(ji,jj) = zp
858               ENDIF
859            END DO
860         END DO
861         CALL lbc_lnk( strength, 'T', 1. )
862         !
863      END SELECT
864      !
865   END SUBROUTINE ice_strength
866
867
868   SUBROUTINE ice_dyn_rdgrft_init
869      !!-------------------------------------------------------------------
870      !!                  ***  ROUTINE ice_dyn_rdgrft_init ***
871      !!
872      !! ** Purpose :   Physical constants and parameters linked
873      !!                to the mechanical ice redistribution
874      !!
875      !! ** Method  :   Read the namdyn_rdgrft namelist
876      !!                and check the parameters values
877      !!                called at the first timestep (nit000)
878      !!
879      !! ** input   :   Namelist namdyn_rdgrft
880      !!-------------------------------------------------------------------
881      INTEGER :: ios                 ! Local integer output status for namelist read
882      !!
883      NAMELIST/namdyn_rdgrft/ ln_str_H79, rn_pstar, rn_crhg, &
884         &                    rn_csrdg  ,                    &
885         &                    ln_partf_lin, rn_gstar,        &
886         &                    ln_partf_exp, rn_astar,        & 
887         &                    ln_ridging, rn_hstar, rn_porordg, rn_fsnwrdg, rn_fpndrdg,  & 
888         &                    ln_rafting, rn_hraft, rn_craft  , rn_fsnwrft, rn_fpndrft
889      !!-------------------------------------------------------------------
890      !
891      REWIND( numnam_ice_ref )              ! Namelist namicetdme in reference namelist : Ice mechanical ice redistribution
892      READ  ( numnam_ice_ref, namdyn_rdgrft, IOSTAT = ios, ERR = 901)
893901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdyn_rdgrft in reference namelist', lwp )
894      REWIND( numnam_ice_cfg )              ! Namelist namdyn_rdgrft in configuration namelist : Ice mechanical ice redistribution
895      READ  ( numnam_ice_cfg, namdyn_rdgrft, IOSTAT = ios, ERR = 902 )
896902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdyn_rdgrft in configuration namelist', lwp )
897      IF(lwm) WRITE ( numoni, namdyn_rdgrft )
898      !
899      IF (lwp) THEN                          ! control print
900         WRITE(numout,*)
901         WRITE(numout,*) 'ice_dyn_rdgrft_init: ice parameters for ridging/rafting '
902         WRITE(numout,*) '~~~~~~~~~~~~~~~~~~'
903         WRITE(numout,*) '   Namelist namdyn_rdgrft:'
904         WRITE(numout,*) '      ice strength parameterization Hibler (1979)              ln_str_H79   = ', ln_str_H79 
905         WRITE(numout,*) '            1st bulk-rheology parameter                        rn_pstar     = ', rn_pstar
906         WRITE(numout,*) '            2nd bulk-rhelogy parameter                         rn_crhg      = ', rn_crhg
907         WRITE(numout,*) '      Fraction of shear energy contributing to ridging         rn_csrdg     = ', rn_csrdg 
908         WRITE(numout,*) '      linear ridging participation function                    ln_partf_lin = ', ln_partf_lin
909         WRITE(numout,*) '            Fraction of ice coverage contributing to ridging   rn_gstar     = ', rn_gstar
910         WRITE(numout,*) '      Exponential ridging participation function               ln_partf_exp = ', ln_partf_exp
911         WRITE(numout,*) '            Equivalent to G* for an exponential function       rn_astar     = ', rn_astar
912         WRITE(numout,*) '      Ridging of ice sheets or not                             ln_ridging   = ', ln_ridging
913         WRITE(numout,*) '            max ridged ice thickness                           rn_hstar     = ', rn_hstar
914         WRITE(numout,*) '            Initial porosity of ridges                         rn_porordg   = ', rn_porordg
915         WRITE(numout,*) '            Fraction of snow volume conserved during ridging   rn_fsnwrdg   = ', rn_fsnwrdg 
916         WRITE(numout,*) '            Fraction of pond volume conserved during ridging   rn_fpndrdg   = ', rn_fpndrdg 
917         WRITE(numout,*) '      Rafting of ice sheets or not                             ln_rafting   = ', ln_rafting
918         WRITE(numout,*) '            Parmeter thickness (threshold between ridge-raft)  rn_hraft     = ', rn_hraft
919         WRITE(numout,*) '            Rafting hyperbolic tangent coefficient             rn_craft     = ', rn_craft 
920         WRITE(numout,*) '            Fraction of snow volume conserved during rafting   rn_fsnwrft   = ', rn_fsnwrft 
921         WRITE(numout,*) '            Fraction of pond volume conserved during rafting   rn_fpndrft   = ', rn_fpndrft 
922      ENDIF
923      !
924      IF ( ( ln_partf_lin .AND. ln_partf_exp ) .OR. ( .NOT.ln_partf_lin .AND. .NOT.ln_partf_exp ) ) THEN
925         CALL ctl_stop( 'ice_dyn_rdgrft_init: choose one and only one participation function (ln_partf_lin or ln_partf_exp)' )
926      ENDIF
927      !                              ! allocate tke arrays
928      IF( ice_dyn_rdgrft_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'ice_dyn_rdgrft_init: unable to allocate arrays' )
929      !
930  END SUBROUTINE ice_dyn_rdgrft_init
931
932#else
933   !!----------------------------------------------------------------------
934   !!   Default option         Empty module           NO ESIM sea-ice model
935   !!----------------------------------------------------------------------
936#endif
937
938   !!======================================================================
939END MODULE icedyn_rdgrft
Note: See TracBrowser for help on using the repository browser.