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

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

change variable names

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