source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerdgrft.F90 @ 8498

Last change on this file since 8498 was 8498, checked in by clem, 3 years ago

changes in style - part2 -

File size: 47.8 KB
Line 
1MODULE icerdgrft
2   !!======================================================================
3   !!                       ***  MODULE icerdgrft ***
4   !! LIM-3 : 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'                                      LIM-3 sea-ice model
13   !!----------------------------------------------------------------------
14   USE par_oce        ! ocean parameters
15   USE dom_oce        ! ocean domain
16   USE phycst         ! physical constants (ocean directory)
17   USE sbc_oce , ONLY : sss_m, sst_m   ! surface boundary condition: ocean fields
18   USE ice1D          ! sea-ice: thermodynamics
19   USE ice            ! sea-ice: variables
20   USE icevar         ! sea-ice: operations
21   USE icectl         ! sea-ice: control prints
22   !
23   USE lbclnk         ! lateral boundary condition - MPP exchanges
24   USE lib_mpp        ! MPP library
25   USE in_out_manager ! I/O manager
26   USE iom            ! I/O manager
27   USE lib_fortran    ! glob_sum
28   USE timing         ! Timing
29
30   IMPLICIT NONE
31   PRIVATE
32
33   PUBLIC   ice_rdgrft               ! called by ice_stp
34   PUBLIC   ice_rdgrft_icestrength   ! called by icerhg_evp
35   PUBLIC   ice_rdgrft_init          ! called by ice_stp
36   PUBLIC   ice_rdgrft_alloc         ! called by ice_init
37
38   !-----------------------------------------------------------------------
39   ! Variables shared among ridging subroutines
40   !-----------------------------------------------------------------------
41   REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   asum     ! sum of total ice and open water area
42   REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   aksum    ! ratio of area removed to area ridged
43   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   athorn   ! participation function; fraction of ridging/closing associated w/ category n
44   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   hrmin    ! minimum ridge thickness
45   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   hrmax    ! maximum ridge thickness
46   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   hraft    ! thickness of rafted ice
47   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   krdg     ! thickness of ridging ice / mean ridge thickness
48   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   aridge   ! participating ice ridging
49   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   araft    ! participating ice rafting
50   !
51   REAL(wp), PARAMETER ::   krdgmin = 1.1_wp    ! min ridge thickness multiplier
52   REAL(wp), PARAMETER ::   kraft   = 0.5_wp    ! rafting multipliyer
53   !
54   REAL(wp)            ::   zdrho               !
55   !
56   !
57   !!----------------------------------------------------------------------
58   !! NEMO/ICE 4.0 , NEMO Consortium (2017)
59   !! $Id: icerdgrft.F90 8378 2017-07-26 13:55:59Z clem $
60   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
61   !!----------------------------------------------------------------------
62CONTAINS
63
64   INTEGER FUNCTION ice_rdgrft_alloc()
65      !!---------------------------------------------------------------------!
66      !!                ***  ROUTINE ice_rdgrft_alloc ***
67      !!---------------------------------------------------------------------!
68      ALLOCATE( asum (jpi,jpj)     , athorn(jpi,jpj,0:jpl) , aksum (jpi,jpj)     ,     &
69         &      hrmin(jpi,jpj,jpl) , hraft (jpi,jpj,jpl)   , aridge(jpi,jpj,jpl) ,     &
70         &      hrmax(jpi,jpj,jpl) , krdg  (jpi,jpj,jpl)   , araft (jpi,jpj,jpl) , STAT=ice_rdgrft_alloc )
71         !
72      IF( ice_rdgrft_alloc /= 0 )   CALL ctl_warn( 'ice_rdgrft_alloc: failed to allocate arrays' )
73      !
74   END FUNCTION ice_rdgrft_alloc
75
76
77   SUBROUTINE ice_rdgrft( kt )
78      !!---------------------------------------------------------------------!
79      !!                ***  ROUTINE ice_rdgrft ***
80      !!
81      !! ** Purpose :   computes the mechanical redistribution of ice thickness
82      !!
83      !! ** Method  :   Steps :
84      !!       1) Thickness categories boundaries, ice / o.w. concentrations
85      !!          Ridge preparation
86      !!       2) Dynamical inputs (closing rate, divu_adv, opning)
87      !!       3) Ridging iteration
88      !!       4) Ridging diagnostics
89      !!       5) Heat, salt and freshwater fluxes
90      !!       6) Compute increments of tate variables and come back to old values
91      !!
92      !! References :   Flato, G. M., and W. D. Hibler III, 1995, JGR, 100, 18,611-18,626.
93      !!                Hibler, W. D. III, 1980, MWR, 108, 1943-1973, 1980.
94      !!                Rothrock, D. A., 1975: JGR, 80, 4514-4519.
95      !!                Thorndike et al., 1975, JGR, 80, 4501-4513.
96      !!                Bitz et al., JGR, 2001
97      !!                Amundrud and Melling, JGR 2005
98      !!                Babko et al., JGR 2002
99      !!
100      !!     This routine is based on CICE code and authors William H. Lipscomb,
101      !!  and Elizabeth C. Hunke, LANL are gratefully acknowledged
102      !!--------------------------------------------------------------------!
103      INTEGER, INTENT(in) ::   kt     ! number of iteration
104      !!
105      INTEGER  ::   ji, jj, jk, jl     ! dummy loop index
106      INTEGER  ::   niter              ! local integer
107      INTEGER  ::   iterate_ridging    ! if =1, repeat the ridging
108      REAL(wp) ::   za                 ! local scalar
109      REAL(wp), DIMENSION(jpi,jpj) ::   closing_net     ! net rate at which area is removed    (1/s)
110      !                                                 ! (ridging ice area - area of new ridges) / dt
111      REAL(wp), DIMENSION(jpi,jpj) ::   divu_adv        ! divu as implied by transport scheme  (1/s)
112      REAL(wp), DIMENSION(jpi,jpj) ::   opning          ! rate of opening due to divergence/shear
113      REAL(wp), DIMENSION(jpi,jpj) ::   closing_gross   ! rate at which area removed, not counting area of new ridges
114      !
115      INTEGER, PARAMETER ::   nitermax = 20   
116      !
117      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
118      !!-----------------------------------------------------------------------------
119      IF( nn_timing == 1 )  CALL timing_start('icerdgrft')
120
121      IF( kt == nit000 ) THEN
122         IF(lwp) WRITE(numout,*)
123         IF(lwp) WRITE(numout,*)'icerdgrft : ice ridging and rafting'
124         IF(lwp) WRITE(numout,*)'~~~~~~~~~~'
125         !
126         zdrho = 0.5 * grav * (rau0-rhoic) * rhoic * r1_rau0      ! proport const for PE
127         !
128      ENDIF     
129      !                    ! conservation test
130      IF( ln_limdiachk )   CALL ice_cons_hsm(0, 'icerdgrft', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
131
132      !-----------------------------------------------------------------------------!
133      ! 1) Thickness categories boundaries, ice / o.w. concentrations, init_ons
134      !-----------------------------------------------------------------------------!
135      !
136      CALL ice_rdgrft_ridgeprep                             ! prepare ridging
137      !
138      DO jj = 1, jpj                                        ! Initialize arrays.
139         DO ji = 1, jpi
140            !-----------------------------------------------------------------------------!
141            ! 2) Dynamical inputs (closing rate, divu_adv, opning)
142            !-----------------------------------------------------------------------------!
143            !
144            ! 2.1 closing_net
145            !-----------------
146            ! Compute the net rate of closing due to convergence
147            ! and shear, based on Flato and Hibler (1995).
148            !
149            ! The energy dissipation rate is equal to the net closing rate
150            ! times the ice strength.
151            !
152            ! NOTE: The NET closing rate is equal to the rate that open water
153            !  area is removed, plus the rate at which ice area is removed by
154            !  ridging, minus the rate at which area is added in new ridges.
155            !  The GROSS closing rate is equal to the first two terms (open
156            !  water closing and thin ice ridging) without the third term
157            !  (thick, newly ridged ice).
158
159            closing_net(ji,jj) = rn_cs * 0.5_wp * ( delta_i(ji,jj) - ABS( divu_i(ji,jj) ) ) - MIN( divu_i(ji,jj), 0._wp )
160
161            ! 2.2 divu_adv
162            !--------------
163            ! Compute divu_adv, the divergence rate given by the transport/
164            ! advection scheme, which may not be equal to divu as computed
165            ! from the velocity field.
166            !
167            ! If divu_adv < 0, make sure the closing rate is large enough
168            ! to give asum = 1.0 after ridging.
169           
170            divu_adv(ji,jj) = ( 1._wp - asum(ji,jj) ) * r1_rdtice  ! asum found in ridgeprep
171
172            IF( divu_adv(ji,jj) < 0._wp )   closing_net(ji,jj) = MAX( closing_net(ji,jj), -divu_adv(ji,jj) )
173
174            ! 2.3 opning
175            !------------
176            ! Compute the (non-negative) opening rate that will give asum = 1.0 after ridging.
177            opning(ji,jj) = closing_net(ji,jj) + divu_adv(ji,jj)
178         END DO
179      END DO
180
181      !-----------------------------------------------------------------------------!
182      ! 3) Ridging iteration
183      !-----------------------------------------------------------------------------!
184      niter           = 1                 ! iteration counter
185      iterate_ridging = 1
186
187      DO WHILE ( iterate_ridging > 0 .AND. niter < nitermax )
188
189         ! 3.2 closing_gross
190         !-----------------------------------------------------------------------------!
191         ! Based on the ITD of ridging and ridged ice, convert the net
192         !  closing rate to a gross closing rate. 
193         ! NOTE: 0 < aksum <= 1
194         closing_gross(:,:) = closing_net(:,:) / aksum(:,:)
195
196         ! correction to closing rate and opening if closing rate is excessive
197         !---------------------------------------------------------------------
198         ! Reduce the closing rate if more than 100% of the open water
199         ! would be removed.  Reduce the opening rate proportionately.
200         DO jj = 1, jpj
201            DO ji = 1, jpi
202               za   = ( opning(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) ) * rdt_ice
203               IF    ( za < 0._wp .AND. za > - ato_i(ji,jj) ) THEN                  ! would lead to negative ato_i
204                  opning(ji,jj) = athorn(ji,jj,0) * closing_gross(ji,jj) - ato_i(ji,jj) * r1_rdtice 
205               ELSEIF( za > 0._wp .AND. za > ( asum(ji,jj) - ato_i(ji,jj) ) ) THEN  ! would lead to ato_i > asum
206                  opning(ji,jj) = athorn(ji,jj,0) * closing_gross(ji,jj) + ( asum(ji,jj) - ato_i(ji,jj) ) * r1_rdtice 
207               ENDIF
208            END DO
209         END DO
210
211         ! correction to closing rate / opening if excessive ice removal
212         !---------------------------------------------------------------
213         ! Reduce the closing rate if more than 100% of any ice category
214         ! would be removed.  Reduce the opening rate proportionately.
215         DO jl = 1, jpl
216            DO jj = 1, jpj
217               DO ji = 1, jpi
218                  za = athorn(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice
219                  IF( za  >  a_i(ji,jj,jl) ) THEN
220                     closing_gross(ji,jj) = closing_gross(ji,jj) * a_i(ji,jj,jl) / za
221                  ENDIF
222               END DO
223            END DO
224         END DO
225
226         ! 3.3 Redistribute area, volume, and energy.
227         !-----------------------------------------------------------------------------!
228         CALL ice_rdgrft_ridgeshift( opning, closing_gross )
229         
230         ! 3.4 Compute total area of ice plus open water after ridging.
231         !-----------------------------------------------------------------------------!
232         ! This is in general not equal to one because of divergence during transport
233         asum(:,:) = ato_i(:,:) + SUM( a_i, dim=3 )
234
235         ! 3.5 Do we keep on iterating?
236         !-----------------------------------------------------------------------------!
237         ! Check whether asum = 1.  If not (because the closing and opening
238         ! rates were reduced above), ridge again with new rates.
239         iterate_ridging = 0
240         DO jj = 1, jpj
241            DO ji = 1, jpi
242               IF( ABS( asum(ji,jj) - 1._wp ) < epsi10 ) THEN
243                  closing_net(ji,jj) = 0._wp
244                  opning     (ji,jj) = 0._wp
245                  ato_i      (ji,jj) = MAX( 0._wp, 1._wp - SUM( a_i(ji,jj,:) ) )
246               ELSE
247                  iterate_ridging    = 1
248                  divu_adv   (ji,jj) = ( 1._wp - asum(ji,jj) ) * r1_rdtice
249                  closing_net(ji,jj) = MAX( 0._wp, -divu_adv(ji,jj) )
250                  opning     (ji,jj) = MAX( 0._wp,  divu_adv(ji,jj) )
251               ENDIF
252            END DO
253         END DO
254         IF( lk_mpp )   CALL mpp_max( iterate_ridging )
255
256         ! Repeat if necessary.
257         ! NOTE: If strength smoothing is turned on, the ridging must be
258         ! iterated globally because of the boundary update in the smoothing.
259         niter = niter + 1
260         !
261         IF( iterate_ridging == 1 ) THEN
262            CALL ice_rdgrft_ridgeprep
263            IF( niter  >  nitermax ) THEN
264               WRITE(numout,*) ' ALERTE : non-converging ridging scheme '
265               WRITE(numout,*) ' niter, iterate_ridging ', niter, iterate_ridging
266            ENDIF
267         ENDIF
268         !
269      END DO !! on the do while over iter
270
271      CALL ice_var_agg( 1 ) 
272
273      !-----------------------------------------------------------------------------!
274      ! control prints
275      !-----------------------------------------------------------------------------!
276      !                    ! conservation test
277      IF( ln_limdiachk )   CALL ice_cons_hsm(1, 'icerdgrft', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
278
279      !                    ! control prints
280      IF( ln_ctl       )   CALL ice_prt3D( 'icerdgrft' )
281      !
282      IF( nn_timing == 1 )  CALL timing_stop('icerdgrft')
283      !
284   END SUBROUTINE ice_rdgrft
285
286
287   SUBROUTINE ice_rdgrft_ridgeprep
288      !!---------------------------------------------------------------------!
289      !!                ***  ROUTINE ice_rdgrft_ridgeprep ***
290      !!
291      !! ** Purpose :   preparation for ridging and strength calculations
292      !!
293      !! ** Method  :   Compute the thickness distribution of the ice and open water
294      !!              participating in ridging and of the resulting ridges.
295      !!---------------------------------------------------------------------!
296      INTEGER  ::   ji, jj, jl    ! dummy loop indices
297      REAL(wp) ::   z1_gstar, z1_astar, zhmean, zdummy   ! local scalar
298      REAL(wp), DIMENSION(jpi,jpj,-1:jpl) ::   zGsum     ! zGsum(n) = sum of areas in categories 0 to n
299      !------------------------------------------------------------------------------!
300
301      z1_gstar      = 1._wp / rn_gstar   
302      z1_astar      = 1._wp / rn_astar   
303
304      CALL ice_var_zapsmall   ! Zero out categories with very small areas
305
306      !                       ! Ice thickness needed for rafting
307      WHERE( a_i(:,:,:) >= epsi20 )   ;   ht_i(:,:,:) = v_i (:,:,:) / a_i(:,:,:)
308      ELSEWHERE                       ;   ht_i(:,:,:) = 0._wp
309      END WHERE
310
311      !------------------------------------------------------------------------------!
312      ! 1) Participation function
313      !------------------------------------------------------------------------------!
314      !
315      ! Compute total area of ice plus open water.
316      ! This is in general not equal to one because of divergence during transport
317      asum(:,:) = ato_i(:,:) + SUM( a_i, dim=3 )
318      !
319      ! Compute cumulative thickness distribution function
320      ! Compute the cumulative thickness distribution function zGsum,
321      ! where zGsum(n) is the fractional area in categories 0 to n.
322      ! initial value (in h = 0) equals open water area
323      zGsum(:,:,-1) = 0._wp
324      zGsum(:,:,0 ) = ato_i(:,:) / asum(:,:)
325      DO jl = 1, jpl
326         zGsum(:,:,jl) = ( ato_i(:,:) + SUM( a_i(:,:,1:jl), dim=3 ) ) / asum(:,:) ! sum(1:jl) is ok (and not jpl)
327      END DO
328
329      ! 1.3 Compute participation function a(h) = b(h).g(h) (athorn)
330      !--------------------------------------------------------------------------------------------------
331      ! Compute the participation function athorn; this is analogous to
332      ! a(h) = b(h)g(h) as defined in Thorndike et al. (1975).
333      ! area lost from category n due to ridging/closing
334      ! athorn(n)   = total area lost due to ridging/closing
335      ! assume b(h) = (2/Gstar) * (1 - G(h)/Gstar).
336      !
337      ! The expressions for athorn are found by integrating b(h)g(h) between
338      ! the category boundaries.
339      ! athorn is always >= 0 and SUM(athorn(0:jpl))=1
340      !-----------------------------------------------------------------
341      !
342      IF( nn_partfun == 0 ) THEN       !--- Linear formulation (Thorndike et al., 1975)
343         DO jl = 0, jpl   
344            DO jj = 1, jpj 
345               DO ji = 1, jpi
346                  IF    ( zGsum(ji,jj,jl)   < rn_gstar ) THEN
347                     athorn(ji,jj,jl) = z1_gstar * ( zGsum(ji,jj,jl) - zGsum(ji,jj,jl-1) ) * &
348                        &                          ( 2._wp - ( zGsum(ji,jj,jl-1) + zGsum(ji,jj,jl) ) * z1_gstar )
349                  ELSEIF( zGsum(ji,jj,jl-1) < rn_gstar ) THEN
350                     athorn(ji,jj,jl) = z1_gstar * ( rn_gstar        - zGsum(ji,jj,jl-1) ) *  &
351                        &                          ( 2._wp - ( zGsum(ji,jj,jl-1) + rn_gstar        ) * z1_gstar )
352                  ELSE
353                     athorn(ji,jj,jl) = 0._wp
354                  ENDIF
355               END DO
356            END DO
357         END DO
358         !
359      ELSE                             !--- Exponential, more stable formulation (Lipscomb et al, 2007)
360         !                       
361         zdummy = 1._wp / ( 1._wp - EXP(-z1_astar) )        ! precompute exponential terms using zGsum as a work array
362         DO jl = -1, jpl
363            zGsum(:,:,jl) = EXP( -zGsum(:,:,jl) * z1_astar ) * zdummy
364         END DO
365         DO jl = 0, jpl
366            athorn(:,:,jl) = zGsum(:,:,jl-1) - zGsum(:,:,jl)
367         END DO
368         !
369      ENDIF
370
371      !                                !--- Ridging and rafting participation concentrations
372      IF( ln_rafting .AND. ln_ridging ) THEN             !- ridging & rafting
373         DO jl = 1, jpl
374            DO jj = 1, jpj 
375               DO ji = 1, jpi
376                  aridge(ji,jj,jl) = ( 1._wp + TANH ( rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) ) ) * 0.5_wp * athorn(ji,jj,jl)
377                  araft (ji,jj,jl) = athorn(ji,jj,jl) - aridge(ji,jj,jl)
378               END DO
379            END DO
380         END DO
381      ELSEIF( ln_ridging .AND. .NOT. ln_rafting ) THEN   !- ridging alone
382         aridge(:,:,:) = athorn(:,:,1:jpl)
383         araft (:,:,:) = 0._wp
384      ELSEIF( ln_rafting .AND. .NOT. ln_ridging ) THEN   !- rafting alone   
385        aridge(:,:,:) = 0._wp
386         araft (:,:,:) = athorn(:,:,1:jpl)
387      ELSE                                               !- no ridging & no rafting
388         aridge(:,:,:) = 0._wp
389         araft (:,:,:) = 0._wp         
390      ENDIF
391
392      !-----------------------------------------------------------------
393      ! 2) Transfer function
394      !-----------------------------------------------------------------
395      ! Compute max and min ridged ice thickness for each ridging category.
396      ! Assume ridged ice is uniformly distributed between hrmin and hrmax.
397      !
398      ! This parameterization is a modified version of Hibler (1980).
399      ! The mean ridging thickness, zhmean, is proportional to hi^(0.5)
400      !  and for very thick ridging ice must be >= krdgmin*hi
401      !
402      ! The minimum ridging thickness, hrmin, is equal to 2*hi
403      !  (i.e., rafting) and for very thick ridging ice is
404      !  constrained by hrmin <= (zhmean + hi)/2.
405      !
406      ! The maximum ridging thickness, hrmax, is determined by
407      !  zhmean and hrmin.
408      !
409      ! These modifications have the effect of reducing the ice strength
410      ! (relative to the Hibler formulation) when very thick ice is
411      ! ridging.
412      !
413      ! aksum = net area removed/ total area removed
414      ! where total area removed = area of ice that ridges
415      !         net area removed = total area removed - area of new ridges
416      !-----------------------------------------------------------------
417
418      aksum(:,:) = athorn(:,:,0)
419      ! Transfer function
420      DO jl = 1, jpl !all categories have a specific transfer function
421         DO jj = 1, jpj
422            DO ji = 1, jpi
423               IF ( athorn(ji,jj,jl) > 0._wp ) THEN
424                  zhmean          = MAX( SQRT( rn_hstar * ht_i(ji,jj,jl) ), ht_i(ji,jj,jl) * krdgmin )
425                  hrmin(ji,jj,jl) = MIN( 2._wp * ht_i(ji,jj,jl), 0.5_wp * ( zhmean + ht_i(ji,jj,jl) ) )
426                  hrmax(ji,jj,jl) = 2._wp * zhmean - hrmin(ji,jj,jl)
427                  hraft(ji,jj,jl) = ht_i(ji,jj,jl) / kraft
428                  krdg (ji,jj,jl) = ht_i(ji,jj,jl) / MAX( zhmean, epsi20 )
429                  !
430                  ! Normalization factor : aksum, ensures mass conservation
431                  aksum(ji,jj) = aksum(ji,jj) + aridge(ji,jj,jl) * ( 1._wp - krdg(ji,jj,jl) )    &
432                     &                        + araft (ji,jj,jl) * ( 1._wp - kraft          )
433               ELSE
434                  hrmin(ji,jj,jl)  = 0._wp 
435                  hrmax(ji,jj,jl)  = 0._wp 
436                  hraft(ji,jj,jl)  = 0._wp 
437                  krdg (ji,jj,jl)  = 1._wp
438               ENDIF
439            END DO
440         END DO
441      END DO
442      !
443   END SUBROUTINE ice_rdgrft_ridgeprep
444
445
446   SUBROUTINE ice_rdgrft_ridgeshift( opning, closing_gross )
447      !!----------------------------------------------------------------------
448      !!                ***  ROUTINE ice_rdgrft_icestrength ***
449      !!
450      !! ** Purpose :   shift ridging ice among thickness categories of ice thickness
451      !!
452      !! ** Method  :   Remove area, volume, and energy from each ridging category
453      !!              and add to thicker ice categories.
454      !!----------------------------------------------------------------------
455      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   opning         ! rate of opening due to divergence/shear
456      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   closing_gross  ! rate at which area retreats, excluding area of new ridges
457      !
458      INTEGER ::   ji, jj, jl, jl1, jl2, jk   ! dummy loop indices
459      INTEGER ::   ij                ! horizontal index, combines i and j loops
460      INTEGER ::   icells            ! number of cells with a_i > puny
461      REAL(wp) ::   hL, hR, farea    ! left and right limits of integration and new area going to jl2
462
463      INTEGER , DIMENSION(jpij) ::   indxi, indxj   ! compressed indices
464      REAL(wp), DIMENSION(jpij) ::   zswitch, fvol   ! new ridge volume going to jl2
465
466      REAL(wp), DIMENSION(jpij) ::   afrac            ! fraction of category area ridged
467      REAL(wp), DIMENSION(jpij) ::   ardg1 , ardg2    ! area of ice ridged & new ridges
468      REAL(wp), DIMENSION(jpij) ::   vsrdg , esrdg    ! snow volume & energy of ridging ice
469      ! MV MP 2016
470      REAL(wp), DIMENSION(jpij) ::   vprdg            ! pond volume of ridging ice
471      REAL(wp), DIMENSION(jpij) ::   aprdg1           ! pond area of ridging ice
472      REAL(wp), DIMENSION(jpij) ::   aprdg2           ! pond area of ridging ice
473      ! END MV MP 2016
474      REAL(wp), DIMENSION(jpij) ::   dhr   , dhr2     ! hrmax - hrmin  &  hrmax^2 - hrmin^2
475
476      REAL(wp), DIMENSION(jpij) ::   vrdg1   ! volume of ice ridged
477      REAL(wp), DIMENSION(jpij) ::   vrdg2   ! volume of new ridges
478      REAL(wp), DIMENSION(jpij) ::   vsw     ! volume of seawater trapped into ridges
479      REAL(wp), DIMENSION(jpij) ::   srdg1   ! sal*volume of ice ridged
480      REAL(wp), DIMENSION(jpij) ::   srdg2   ! sal*volume of new ridges
481      REAL(wp), DIMENSION(jpij) ::   smsw    ! sal*volume of water trapped into ridges
482      REAL(wp), DIMENSION(jpij) ::   oirdg1, oirdg2   ! ice age of ice ridged
483
484      REAL(wp), DIMENSION(jpij) ::   afrft            ! fraction of category area rafted
485      REAL(wp), DIMENSION(jpij) ::   arft1 , arft2    ! area of ice rafted and new rafted zone
486      REAL(wp), DIMENSION(jpij) ::   virft , vsrft    ! ice & snow volume of rafting ice
487      ! MV MP 2016
488      REAL(wp), DIMENSION(jpij) ::   vprft            ! pond volume of rafting ice
489      REAL(wp), DIMENSION(jpij) ::   aprft1           ! pond area of rafted ice
490      REAL(wp), DIMENSION(jpij) ::   aprft2           ! pond area of new rafted ice
491      ! END MV MP 2016
492      REAL(wp), DIMENSION(jpij) ::   esrft , smrft    ! snow energy & salinity of rafting ice
493      REAL(wp), DIMENSION(jpij) ::   oirft1, oirft2   ! ice age of ice rafted
494
495      REAL(wp), DIMENSION(jpij,nlay_i) ::   eirft      ! ice energy of rafting ice
496      REAL(wp), DIMENSION(jpij,nlay_i) ::   erdg1      ! enth*volume of ice ridged
497      REAL(wp), DIMENSION(jpij,nlay_i) ::   erdg2      ! enth*volume of new ridges
498      REAL(wp), DIMENSION(jpij,nlay_i) ::   ersw       ! enth of water trapped into ridges
499      !!----------------------------------------------------------------------
500
501      !-------------------------------------------------------------------------------
502      ! 1) Compute change in open water area due to closing and opening.
503      !-------------------------------------------------------------------------------
504      DO jj = 1, jpj
505         DO ji = 1, jpi
506            ato_i(ji,jj) = MAX( 0._wp, ato_i(ji,jj) +  &
507               &                     ( opning(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) ) * rdt_ice )
508         END DO
509      END DO
510
511      !-----------------------------------------------------------------
512      ! 3) Pump everything from ice which is being ridged / rafted
513      !-----------------------------------------------------------------
514      ! Compute the area, volume, and energy of ice ridging in each
515      ! category, along with the area of the resulting ridge.
516
517      DO jl1 = 1, jpl !jl1 describes the ridging category
518
519         !------------------------------------------------
520         ! 3.1) Identify grid cells with nonzero ridging
521         !------------------------------------------------
522         icells = 0
523         DO jj = 1, jpj
524            DO ji = 1, jpi
525               IF( athorn(ji,jj,jl1) > 0._wp .AND. closing_gross(ji,jj) > 0._wp ) THEN
526                  icells = icells + 1
527                  indxi(icells) = ji
528                  indxj(icells) = jj
529               ENDIF
530            END DO
531         END DO
532
533         DO ij = 1, icells
534            ji = indxi(ij) ; jj = indxj(ij)
535
536            !--------------------------------------------------------------------
537            ! 3.2) Compute area of ridging ice (ardg1) and of new ridge (ardg2)
538            !--------------------------------------------------------------------
539            ardg1(ij) = aridge(ji,jj,jl1) * closing_gross(ji,jj) * rdt_ice
540            arft1(ij) = araft (ji,jj,jl1) * closing_gross(ji,jj) * rdt_ice
541
542            !---------------------------------------------------------------
543            ! 3.3) Compute ridging /rafting fractions, make sure afrac <=1
544            !---------------------------------------------------------------
545            afrac(ij) = ardg1(ij) / a_i(ji,jj,jl1) !ridging
546            afrft(ij) = arft1(ij) / a_i(ji,jj,jl1) !rafting
547            ardg2(ij) = ardg1(ij) * krdg(ji,jj,jl1)
548            arft2(ij) = arft1(ij) * kraft
549
550            !--------------------------------------------------------------------------
551            ! 3.4) Substract area, volume, and energy from ridging
552            !     / rafting category n1.
553            !--------------------------------------------------------------------------
554            vrdg1(ij) = v_i(ji,jj,jl1) * afrac(ij)
555            vrdg2(ij) = vrdg1(ij) * ( 1. + rn_por_rdg )
556            vsw  (ij) = vrdg1(ij) * rn_por_rdg
557
558            vsrdg (ij) = v_s  (ji,jj,  jl1) * afrac(ij)
559            esrdg (ij) = e_s  (ji,jj,1,jl1) * afrac(ij)
560            !MV MP 2016
561            IF ( nn_pnd_scheme > 0 ) THEN
562               aprdg1(ij) = a_ip(ji,jj, jl1) * afrac(ij)
563               aprdg2(ij) = a_ip(ji,jj, jl1) * afrac(ij) * krdg(ji,jj,jl1)
564               vprdg(ij)  = v_ip(ji,jj, jl1) * afrac(ij)
565            ENDIF
566            ! END MV MP 2016
567            srdg1 (ij) = smv_i(ji,jj,  jl1) * afrac(ij)
568            oirdg1(ij) = oa_i (ji,jj,  jl1) * afrac(ij)
569            oirdg2(ij) = oa_i (ji,jj,  jl1) * afrac(ij) * krdg(ji,jj,jl1) 
570
571            ! rafting volumes, heat contents ...
572            virft (ij) = v_i  (ji,jj,  jl1) * afrft(ij)
573            vsrft (ij) = v_s  (ji,jj,  jl1) * afrft(ij)
574            !MV MP 2016
575            IF ( nn_pnd_scheme > 0 ) THEN
576               aprft1(ij) = a_ip (ji,jj,  jl1) * afrft(ij)
577               aprft2(ij) = a_ip (ji,jj,  jl1) * afrft(ij) * kraft
578               vprft(ij)  = v_ip(ji,jj,jl1)    * afrft(ij)
579            ENDIF
580            ! END MV MP 2016
581            srdg1 (ij) = smv_i(ji,jj,  jl1) * afrac(ij)
582            esrft (ij) = e_s  (ji,jj,1,jl1) * afrft(ij)
583            smrft (ij) = smv_i(ji,jj,  jl1) * afrft(ij) 
584            oirft1(ij) = oa_i (ji,jj,  jl1) * afrft(ij) 
585            oirft2(ij) = oa_i (ji,jj,  jl1) * afrft(ij) * kraft 
586
587            !-----------------------------------------------------------------
588            ! 3.5) Compute properties of new ridges
589            !-----------------------------------------------------------------
590            smsw(ij)  = vsw(ij) * sss_m(ji,jj)                   ! salt content of seawater frozen in voids
591            srdg2(ij) = srdg1(ij) + smsw(ij)                     ! salt content of new ridge
592           
593            sfx_dyn(ji,jj) = sfx_dyn(ji,jj) - smsw(ij) * rhoic * r1_rdtice
594            wfx_dyn(ji,jj) = wfx_dyn(ji,jj) - vsw (ij) * rhoic * r1_rdtice   ! increase in ice volume due to seawater frozen in voids
595
596            ! virtual salt flux to keep salinity constant
597            IF( nn_icesal == 1 .OR. nn_icesal == 3 )  THEN
598               srdg2(ij)      = srdg2(ij) - vsw(ij) * ( sss_m(ji,jj) - sm_i(ji,jj,jl1) )           ! ridge salinity = sm_i
599               sfx_bri(ji,jj) = sfx_bri(ji,jj) + sss_m(ji,jj)    * vsw(ij) * rhoic * r1_rdtice  &  ! put back sss_m into the ocean
600                  &                            - sm_i(ji,jj,jl1) * vsw(ij) * rhoic * r1_rdtice     ! and get  sm_i  from the ocean
601            ENDIF
602
603            !------------------------------------------           
604            ! 3.7 Put the snow somewhere in the ocean
605            !------------------------------------------           
606            !  Place part of the snow lost by ridging into the ocean.
607            !  Note that esrdg > 0; the ocean must cool to melt snow.
608            !  If the ocean temp = Tf already, new ice must grow.
609            !  During the next time step, thermo_rates will determine whether
610            !  the ocean cools or new ice grows.
611            wfx_snw_dyn(ji,jj) = wfx_snw_dyn(ji,jj) + ( rhosn * vsrdg(ij) * ( 1._wp - rn_fsnowrdg )   & 
612               &                                      + rhosn * vsrft(ij) * ( 1._wp - rn_fsnowrft ) ) * r1_rdtice  ! fresh water source for ocean
613
614!!!            wfx_snw_dyn(ji,jj) = wfx_snw_dyn(ji,jj) + zwfx_snw
615
616            hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ( - esrdg(ij) * ( 1._wp - rn_fsnowrdg )         & 
617               &                                - esrft(ij) * ( 1._wp - rn_fsnowrft ) ) * r1_rdtice        ! heat sink for ocean (<0, W.m-2)
618
619            ! MV MP 2016
620            !------------------------------------------           
621            ! 3.X Put the melt pond water in the ocean
622            !------------------------------------------           
623            !  Place part of the melt pond volume into the ocean.
624            IF ( ( nn_pnd_scheme > 0 ) .AND. ln_pnd_fw ) THEN
625               wfx_pnd(ji,jj) = wfx_pnd(ji,jj) + ( rhofw * vprdg(ij) * ( 1._wp - rn_fpondrdg )   & 
626               &                                 + rhofw * vprft(ij) * ( 1._wp - rn_fpondrft ) ) * r1_rdtice  ! fresh water source for ocean
627            ENDIF
628            ! END MV MP 2016
629
630            !-----------------------------------------------------------------
631            ! 3.8 Compute quantities used to apportion ice among categories
632            ! in the n2 loop below
633            !-----------------------------------------------------------------
634            dhr (ij) = 1._wp / ( hrmax(ji,jj,jl1)                    - hrmin(ji,jj,jl1)                    )
635            dhr2(ij) = 1._wp / ( hrmax(ji,jj,jl1) * hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) * hrmin(ji,jj,jl1) )
636
637
638            ! update jl1 (removing ridged/rafted area)
639            a_i  (ji,jj,  jl1) = a_i  (ji,jj,  jl1) - ardg1 (ij) - arft1 (ij)
640            v_i  (ji,jj,  jl1) = v_i  (ji,jj,  jl1) - vrdg1 (ij) - virft (ij)
641            v_s  (ji,jj,  jl1) = v_s  (ji,jj,  jl1) - vsrdg (ij) - vsrft (ij)
642            e_s  (ji,jj,1,jl1) = e_s  (ji,jj,1,jl1) - esrdg (ij) - esrft (ij)
643            smv_i(ji,jj,  jl1) = smv_i(ji,jj,  jl1) - srdg1 (ij) - smrft (ij)
644            oa_i (ji,jj,  jl1) = oa_i (ji,jj,  jl1) - oirdg1(ij) - oirft1(ij)
645
646            ! MV MP 2016
647            IF ( nn_pnd_scheme > 0 ) THEN
648               v_ip (ji,jj,jl1) = v_ip (ji,jj,jl1) - vprdg (ij) - vprft (ij)
649               a_ip (ji,jj,jl1) = a_ip (ji,jj,jl1) - aprdg1(ij) - aprft1(ij)
650            ENDIF
651            ! END MV MP 2016
652
653         END DO
654
655         !--------------------------------------------------------------------
656         ! 3.9 Compute ridging ice enthalpy, remove it from ridging ice and
657         !      compute ridged ice enthalpy
658         !--------------------------------------------------------------------
659         DO jk = 1, nlay_i
660            DO ij = 1, icells
661               ji = indxi(ij) ; jj = indxj(ij)
662               ! heat content of ridged ice
663               erdg1(ij,jk) = e_i(ji,jj,jk,jl1) * afrac(ij) 
664               eirft(ij,jk) = e_i(ji,jj,jk,jl1) * afrft(ij)               
665               
666               ! enthalpy of the trapped seawater (J/m2, >0)
667               ! clem: if sst>0, then ersw <0 (is that possible?)
668               ersw(ij,jk)  = - rhoic * vsw(ij) * rcp * sst_m(ji,jj) * r1_nlay_i
669
670               ! heat flux to the ocean
671               hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ersw(ij,jk) * r1_rdtice  ! > 0 [W.m-2] ocean->ice flux
672
673               ! it is added to sea ice because the sign convention is the opposite of the sign convention for the ocean
674               erdg2(ij,jk) = erdg1(ij,jk) + ersw(ij,jk)
675
676               ! update jl1
677               e_i(ji,jj,jk,jl1) = e_i(ji,jj,jk,jl1) - erdg1(ij,jk) - eirft(ij,jk)
678
679            END DO
680         END DO
681
682         !-------------------------------------------------------------------------------
683         ! 4) Add area, volume, and energy of new ridge to each category jl2
684         !-------------------------------------------------------------------------------
685         DO jl2  = 1, jpl 
686            ! over categories to which ridged/rafted ice is transferred
687            DO ij = 1, icells
688               ji = indxi(ij) ; jj = indxj(ij)
689
690               ! Compute the fraction of ridged ice area and volume going to thickness category jl2.
691               IF( hrmin(ji,jj,jl1) <= hi_max(jl2) .AND. hrmax(ji,jj,jl1) > hi_max(jl2-1) ) THEN
692                  hL = MAX( hrmin(ji,jj,jl1), hi_max(jl2-1) )
693                  hR = MIN( hrmax(ji,jj,jl1), hi_max(jl2)   )
694                  farea    = ( hR      - hL      ) * dhr(ij) 
695                  fvol(ij) = ( hR * hR - hL * hL ) * dhr2(ij)
696               ELSE
697                  farea    = 0._wp 
698                  fvol(ij) = 0._wp                 
699               ENDIF
700
701               ! Compute the fraction of rafted ice area and volume going to thickness category jl2
702!!gm see above               IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. hraft(ji,jj,jl1) >  hi_max(jl2-1) ) THEN
703               IF( hi_max(jl2-1) < hraft(ji,jj,jl1) .AND. hraft(ji,jj,jl1) <= hi_max(jl2)  ) THEN   ;   zswitch(ij) = 1._wp
704               ELSE                                                                                 ;   zswitch(ij) = 0._wp                 
705               ENDIF
706               !
707               a_i  (ji,jj  ,jl2) = a_i  (ji,jj  ,jl2) + ( ardg2 (ij) * farea    + arft2 (ij) * zswitch(ij) )
708               oa_i (ji,jj  ,jl2) = oa_i (ji,jj  ,jl2) + ( oirdg2(ij) * farea    + oirft2(ij) * zswitch(ij) )
709               v_i  (ji,jj  ,jl2) = v_i  (ji,jj  ,jl2) + ( vrdg2 (ij) * fvol(ij) + virft (ij) * zswitch(ij) )
710               smv_i(ji,jj  ,jl2) = smv_i(ji,jj  ,jl2) + ( srdg2 (ij) * fvol(ij) + smrft (ij) * zswitch(ij) )
711               v_s  (ji,jj  ,jl2) = v_s  (ji,jj  ,jl2) + ( vsrdg (ij) * rn_fsnowrdg * fvol(ij)  +  &
712                  &                                        vsrft (ij) * rn_fsnowrft * zswitch(ij) )
713               e_s  (ji,jj,1,jl2) = e_s  (ji,jj,1,jl2) + ( esrdg (ij) * rn_fsnowrdg * fvol(ij)  +  &
714                  &                                        esrft (ij) * rn_fsnowrft * zswitch(ij) )
715               ! MV MP 2016
716               IF ( nn_pnd_scheme > 0 ) THEN
717                  v_ip (ji,jj,jl2) = v_ip(ji,jj,jl2) + (   vprdg (ij) * rn_fpondrdg * fvol   (ij)   &
718                     &                                   + vprft (ij) * rn_fpondrft * zswitch(ij)   )
719                  a_ip (ji,jj,jl2) = a_ip(ji,jj,jl2) + (   aprdg2(ij) * rn_fpondrdg * farea         & 
720                     &                                   + aprft2(ij) * rn_fpondrft * zswitch(ji)   )
721               ENDIF
722               ! END MV MP 2016
723            END DO
724
725            ! Transfer ice energy to category jl2 by ridging
726            DO jk = 1, nlay_i
727               DO ij = 1, icells
728                  ji = indxi(ij) ; jj = indxj(ij)
729                  e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + erdg2(ij,jk) * fvol(ij) + eirft(ij,jk) * zswitch(ij)                 
730               END DO
731            END DO
732            !
733         END DO ! jl2
734         !
735      END DO ! jl1 (deforming categories)
736
737      ! SIMIP diagnostics
738      diag_dmi_dyn(:,:) = - wfx_dyn(:,:)     + rhoic * diag_trp_vi(:,:)
739      diag_dms_dyn(:,:) = - wfx_snw_dyn(:,:) + rhosn * diag_trp_vs(:,:)     
740      !
741   END SUBROUTINE ice_rdgrft_ridgeshift
742
743
744   SUBROUTINE ice_rdgrft_icestrength( kstrngth )
745      !!----------------------------------------------------------------------
746      !!                ***  ROUTINE ice_rdgrft_icestrength ***
747      !!
748      !! ** Purpose :   computes ice strength used in dynamics routines of ice thickness
749      !!
750      !! ** Method  :   Compute the strength of the ice pack, defined as the energy (J m-2)
751      !!              dissipated per unit area removed from the ice pack under compression,
752      !!              and assumed proportional to the change in potential energy caused
753      !!              by ridging. Note that only Hibler's formulation is stable and that
754      !!              ice strength has to be smoothed
755      !!
756      !! ** Inputs / Ouputs : kstrngth (what kind of ice strength we are using)
757      !!----------------------------------------------------------------------
758      INTEGER, INTENT(in) ::   kstrngth    ! = 1 for Rothrock formulation, 0 for Hibler (1979)
759      !
760      INTEGER             ::   ji,jj, jl   ! dummy loop indices
761      INTEGER             ::   ismooth     ! smoothing the resistance to deformation
762      INTEGER             ::   itframe     ! number of time steps for the P smoothing
763      REAL(wp)            ::   zp, z1_3    ! local scalars
764      REAL(wp), DIMENSION(jpi,jpj) ::   zworka           ! temporary array used here
765      REAL(wp), DIMENSION(jpi,jpj) ::   zstrp1, zstrp2   ! strength at previous time steps
766      !!----------------------------------------------------------------------
767
768      !                              !--------------------------------------------------!
769      CALL ice_rdgrft_ridgeprep      ! Thickness distribution of ridging and ridged ice !
770      !                              !--------------------------------------------------!
771
772      !                              !--------------------------------------------------!
773      IF( kstrngth == 1 ) THEN       ! Ice strength => Rothrock (1975) method           !
774      !                              !--------------------------------------------------!
775         z1_3 = 1._wp / 3._wp
776         DO jl = 1, jpl
777            WHERE( athorn(:,:,jl) > 0._wp )
778               strength(:,:) =  -         athorn(:,:,jl) * ht_i(:,:,jl) * ht_i(:,:,jl)   &  ! PE loss from deforming ice
779                  &             + 2._wp * araft (:,:,jl) * ht_i(:,:,jl) * ht_i(:,:,jl)   &  ! PE gain from rafting ice
780                  &             +         aridge(:,:,jl) * krdg(:,:,jl) * z1_3 *   &        ! PE gain from ridging ice
781                  &                      ( hrmax(:,:,jl) * hrmax(:,:,jl) +         &
782                  &                        hrmin(:,:,jl) * hrmin(:,:,jl) +         &
783                  &                        hrmax(:,:,jl) * hrmin(:,:,jl) )
784            ELSEWHERE
785               strength(:,:) = 0._wp
786            END WHERE
787         END DO
788         strength(:,:) = rn_pe_rdg * zdrho * strength(:,:) / aksum(:,:) * tmask(:,:,1)
789                         ! where zdrho = (g/2)*(rhow-rhoi)*(rhoi/rhow) and rn_pe_rdg accounts for frictional dissipation
790         ismooth = 1
791      !                              !--------------------------------------------------!
792      ELSE                           ! Ice strength => Hibler (1979) method             !
793      !                              !--------------------------------------------------!
794         strength(:,:) = rn_pstar * vt_i(:,:) * EXP( - rn_crhg * ( 1._wp - at_i(:,:) )  ) * tmask(:,:,1)
795         !
796         ismooth = 1
797         !
798      ENDIF
799      !                              !--------------------------------------------------!
800      IF( ln_icestr_bvf ) THEN       ! Impact of brine volume                           !
801      !                              !--------------------------------------------------!
802      ! CAN BE REMOVED
803         DO jj = 1, jpj
804            DO ji = 1, jpi
805               strength(ji,jj) = strength(ji,jj) * exp(-5.88*SQRT(MAX(bvm_i(ji,jj),0.0)))
806            END DO
807         END DO
808      ENDIF
809      !                              !--------------------------------------------------!
810      SELECT CASE( ismooth )         ! Smoothing ice strength                           !
811      !                              !--------------------------------------------------!
812      !
813      CASE( 1 )               !--- Spatial smoothing
814         DO jj = 2, jpjm1
815            DO ji = 2, jpim1
816               IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp ) THEN
817                  zworka(ji,jj) = ( 4.0 * strength(ji,jj)              &
818                     &                  + strength(ji-1,jj) * tmask(ji-1,jj,1) + strength(ji+1,jj) * tmask(ji+1,jj,1) & 
819                     &                  + strength(ji,jj-1) * tmask(ji,jj-1,1) + strength(ji,jj+1) * tmask(ji,jj+1,1) &
820                     &            ) / ( 4.0 + tmask(ji-1,jj,1) + tmask(ji+1,jj,1) + tmask(ji,jj-1,1) + tmask(ji,jj+1,1) )
821               ELSE
822                  zworka(ji,jj) = 0._wp
823               ENDIF
824            END DO
825         END DO
826         
827         DO jj = 2, jpjm1
828            DO ji = 2, jpim1
829               strength(ji,jj) = zworka(ji,jj)
830            END DO
831         END DO
832         CALL lbc_lnk( strength, 'T', 1. )
833         !
834      CASE( 2 )               !--- Temporal smoothing
835         IF ( kt_ice == nit000 ) THEN
836            zstrp1(:,:) = 0._wp
837            zstrp2(:,:) = 0._wp
838         ENDIF
839         !
840         DO jj = 2, jpjm1
841            DO ji = 2, jpim1
842               IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp ) THEN
843                  itframe = 1 ! number of time steps for the running mean
844                  IF ( zstrp1(ji,jj) > 0._wp ) itframe = itframe + 1
845                  IF ( zstrp2(ji,jj) > 0._wp ) itframe = itframe + 1
846                  zp = ( strength(ji,jj) + zstrp1(ji,jj) + zstrp2(ji,jj) ) / itframe
847                  zstrp2  (ji,jj) = zstrp1  (ji,jj)
848                  zstrp1  (ji,jj) = strength(ji,jj)
849                  strength(ji,jj) = zp
850               ENDIF
851            END DO
852         END DO
853         CALL lbc_lnk( strength, 'T', 1. )
854         !
855      END SELECT
856      !
857   END SUBROUTINE ice_rdgrft_icestrength
858
859
860   SUBROUTINE ice_rdgrft_init
861      !!-------------------------------------------------------------------
862      !!                  ***  ROUTINE ice_rdgrft_init ***
863      !!
864      !! ** Purpose :   Physical constants and parameters linked
865      !!                to the mechanical ice redistribution
866      !!
867      !! ** Method  :   Read the namiceitdme namelist
868      !!                and check the parameters values
869      !!                called at the first timestep (nit000)
870      !!
871      !! ** input   :   Namelist namiceitdme
872      !!-------------------------------------------------------------------
873      INTEGER :: ios                 ! Local integer output status for namelist read
874      !!
875      NAMELIST/namiceitdme/ rn_cs     , nn_partfun, rn_gstar  , rn_astar   ,                & 
876         &                  ln_ridging, rn_hstar  , rn_por_rdg, rn_fsnowrdg, rn_fpondrdg,   & 
877         &                  ln_rafting, rn_hraft  , rn_craft  , rn_fsnowrft, rn_fpondrft
878      !!-------------------------------------------------------------------
879      !
880      REWIND( numnam_ice_ref )              ! Namelist namicetdme in reference namelist : Ice mechanical ice redistribution
881      READ  ( numnam_ice_ref, namiceitdme, IOSTAT = ios, ERR = 901)
882901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitdme in reference namelist', lwp )
883      !
884      REWIND( numnam_ice_cfg )              ! Namelist namiceitdme in configuration namelist : Ice mechanical ice redistribution
885      READ  ( numnam_ice_cfg, namiceitdme, IOSTAT = ios, ERR = 902 )
886902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitdme in configuration namelist', lwp )
887      IF(lwm) WRITE ( numoni, namiceitdme )
888      !
889      IF (lwp) THEN                          ! control print
890         WRITE(numout,*)
891         WRITE(numout,*) 'ice_rdgrft_init : ice parameters for mechanical ice redistribution '
892         WRITE(numout,*) '~~~~~~~~~~~~~~~'
893         WRITE(numout,*) '   Namelist namiceitdme'
894         WRITE(numout,*) '      Fraction of shear energy contributing to ridging        rn_cs       = ', rn_cs 
895         WRITE(numout,*) '      Switch for part. function (0) linear (1) exponential    nn_partfun  = ', nn_partfun
896         WRITE(numout,*) '      Fraction of total ice coverage contributing to ridging  rn_gstar    = ', rn_gstar
897         WRITE(numout,*) '      Equivalent to G* for an exponential part function       rn_astar    = ', rn_astar
898         WRITE(numout,*) '      Ridging of ice sheets or not                            ln_ridging  = ', ln_ridging
899         WRITE(numout,*) '      Quantity playing a role in max ridged ice thickness     rn_hstar    = ', rn_hstar
900         WRITE(numout,*) '      Initial porosity of ridges                              rn_por_rdg  = ', rn_por_rdg
901         WRITE(numout,*) '      Fraction of snow volume conserved during ridging        rn_fsnowrdg = ', rn_fsnowrdg 
902         WRITE(numout,*) '      Fraction of pond volume conserved during ridging        rn_fpondrdg = ', rn_fpondrdg 
903         WRITE(numout,*) '      Rafting of ice sheets or not                            ln_rafting  = ', ln_rafting
904         WRITE(numout,*) '      Parmeter thickness (threshold between ridge-raft)       rn_hraft    = ', rn_hraft
905         WRITE(numout,*) '      Rafting hyperbolic tangent coefficient                  rn_craft    = ', rn_craft 
906         WRITE(numout,*) '      Fraction of snow volume conserved during ridging        rn_fsnowrft = ', rn_fsnowrft 
907         WRITE(numout,*) '      Fraction of pond volume conserved during rafting        rn_fpondrft = ', rn_fpondrft 
908      ENDIF
909      !
910   END SUBROUTINE ice_rdgrft_init
911
912#else
913   !!----------------------------------------------------------------------
914   !!   Default option         Empty module          NO LIM-3 sea-ice model
915   !!----------------------------------------------------------------------
916#endif
917
918   !!======================================================================
919END MODULE icerdgrft
Note: See TracBrowser for help on using the repository browser.