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

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

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