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

source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rdgrft.F90 @ 8555

Last change on this file since 8555 was 8555, checked in by clem, 7 years ago

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