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

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

changes in style - part3 - move output into proper files and correct a bug in debug mode

File size: 47.5 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            hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ( - esrdg(ij) * ( 1._wp - rn_fsnowrdg )         & 
615               &                                - esrft(ij) * ( 1._wp - rn_fsnowrft ) ) * r1_rdtice        ! heat sink for ocean (<0, W.m-2)
616
617            ! MV MP 2016
618            !------------------------------------------           
619            ! 3.X Put the melt pond water in the ocean
620            !------------------------------------------           
621            !  Place part of the melt pond volume into the ocean.
622            IF ( ( nn_pnd_scheme > 0 ) .AND. ln_pnd_fw ) THEN
623               wfx_pnd(ji,jj) = wfx_pnd(ji,jj) + ( rhofw * vprdg(ij) * ( 1._wp - rn_fpondrdg )   & 
624               &                                 + rhofw * vprft(ij) * ( 1._wp - rn_fpondrft ) ) * r1_rdtice  ! fresh water source for ocean
625            ENDIF
626            ! END MV MP 2016
627
628            !-----------------------------------------------------------------
629            ! 3.8 Compute quantities used to apportion ice among categories
630            ! in the n2 loop below
631            !-----------------------------------------------------------------
632            dhr (ij) = 1._wp / ( hrmax(ji,jj,jl1)                    - hrmin(ji,jj,jl1)                    )
633            dhr2(ij) = 1._wp / ( hrmax(ji,jj,jl1) * hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) * hrmin(ji,jj,jl1) )
634
635
636            ! update jl1 (removing ridged/rafted area)
637            a_i  (ji,jj,  jl1) = a_i  (ji,jj,  jl1) - ardg1 (ij) - arft1 (ij)
638            v_i  (ji,jj,  jl1) = v_i  (ji,jj,  jl1) - vrdg1 (ij) - virft (ij)
639            v_s  (ji,jj,  jl1) = v_s  (ji,jj,  jl1) - vsrdg (ij) - vsrft (ij)
640            e_s  (ji,jj,1,jl1) = e_s  (ji,jj,1,jl1) - esrdg (ij) - esrft (ij)
641            smv_i(ji,jj,  jl1) = smv_i(ji,jj,  jl1) - srdg1 (ij) - smrft (ij)
642            oa_i (ji,jj,  jl1) = oa_i (ji,jj,  jl1) - oirdg1(ij) - oirft1(ij)
643
644            ! MV MP 2016
645            IF ( nn_pnd_scheme > 0 ) THEN
646               v_ip (ji,jj,jl1) = v_ip (ji,jj,jl1) - vprdg (ij) - vprft (ij)
647               a_ip (ji,jj,jl1) = a_ip (ji,jj,jl1) - aprdg1(ij) - aprft1(ij)
648            ENDIF
649            ! END MV MP 2016
650
651         END DO
652
653         !--------------------------------------------------------------------
654         ! 3.9 Compute ridging ice enthalpy, remove it from ridging ice and
655         !      compute ridged ice enthalpy
656         !--------------------------------------------------------------------
657         DO jk = 1, nlay_i
658            DO ij = 1, icells
659               ji = indxi(ij) ; jj = indxj(ij)
660               ! heat content of ridged ice
661               erdg1(ij,jk) = e_i(ji,jj,jk,jl1) * afrac(ij) 
662               eirft(ij,jk) = e_i(ji,jj,jk,jl1) * afrft(ij)               
663               
664               ! enthalpy of the trapped seawater (J/m2, >0)
665               ! clem: if sst>0, then ersw <0 (is that possible?)
666               ersw(ij,jk)  = - rhoic * vsw(ij) * rcp * sst_m(ji,jj) * r1_nlay_i
667
668               ! heat flux to the ocean
669               hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ersw(ij,jk) * r1_rdtice  ! > 0 [W.m-2] ocean->ice flux
670
671               ! it is added to sea ice because the sign convention is the opposite of the sign convention for the ocean
672               erdg2(ij,jk) = erdg1(ij,jk) + ersw(ij,jk)
673
674               ! update jl1
675               e_i(ji,jj,jk,jl1) = e_i(ji,jj,jk,jl1) - erdg1(ij,jk) - eirft(ij,jk)
676
677            END DO
678         END DO
679
680         !-------------------------------------------------------------------------------
681         ! 4) Add area, volume, and energy of new ridge to each category jl2
682         !-------------------------------------------------------------------------------
683         DO jl2  = 1, jpl 
684            ! over categories to which ridged/rafted ice is transferred
685            DO ij = 1, icells
686               ji = indxi(ij) ; jj = indxj(ij)
687
688               ! Compute the fraction of ridged ice area and volume going to thickness category jl2.
689               IF( hrmin(ji,jj,jl1) <= hi_max(jl2) .AND. hrmax(ji,jj,jl1) > hi_max(jl2-1) ) THEN
690                  hL = MAX( hrmin(ji,jj,jl1), hi_max(jl2-1) )
691                  hR = MIN( hrmax(ji,jj,jl1), hi_max(jl2)   )
692                  farea    = ( hR      - hL      ) * dhr(ij) 
693                  fvol(ij) = ( hR * hR - hL * hL ) * dhr2(ij)
694               ELSE
695                  farea    = 0._wp 
696                  fvol(ij) = 0._wp                 
697               ENDIF
698
699               ! Compute the fraction of rafted ice area and volume going to thickness category jl2
700!!gm see above               IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. hraft(ji,jj,jl1) >  hi_max(jl2-1) ) THEN
701               IF( hi_max(jl2-1) < hraft(ji,jj,jl1) .AND. hraft(ji,jj,jl1) <= hi_max(jl2)  ) THEN   ;   zswitch(ij) = 1._wp
702               ELSE                                                                                 ;   zswitch(ij) = 0._wp                 
703               ENDIF
704               !
705               a_i  (ji,jj  ,jl2) = a_i  (ji,jj  ,jl2) + ( ardg2 (ij) * farea    + arft2 (ij) * zswitch(ij) )
706               oa_i (ji,jj  ,jl2) = oa_i (ji,jj  ,jl2) + ( oirdg2(ij) * farea    + oirft2(ij) * zswitch(ij) )
707               v_i  (ji,jj  ,jl2) = v_i  (ji,jj  ,jl2) + ( vrdg2 (ij) * fvol(ij) + virft (ij) * zswitch(ij) )
708               smv_i(ji,jj  ,jl2) = smv_i(ji,jj  ,jl2) + ( srdg2 (ij) * fvol(ij) + smrft (ij) * zswitch(ij) )
709               v_s  (ji,jj  ,jl2) = v_s  (ji,jj  ,jl2) + ( vsrdg (ij) * rn_fsnowrdg * fvol(ij)  +  &
710                  &                                        vsrft (ij) * rn_fsnowrft * zswitch(ij) )
711               e_s  (ji,jj,1,jl2) = e_s  (ji,jj,1,jl2) + ( esrdg (ij) * rn_fsnowrdg * fvol(ij)  +  &
712                  &                                        esrft (ij) * rn_fsnowrft * zswitch(ij) )
713               ! MV MP 2016
714               IF ( nn_pnd_scheme > 0 ) THEN
715                  v_ip (ji,jj,jl2) = v_ip(ji,jj,jl2) + (   vprdg (ij) * rn_fpondrdg * fvol   (ij)   &
716                     &                                   + vprft (ij) * rn_fpondrft * zswitch(ij)   )
717                  a_ip (ji,jj,jl2) = a_ip(ji,jj,jl2) + (   aprdg2(ij) * rn_fpondrdg * farea         & 
718                     &                                   + aprft2(ij) * rn_fpondrft * zswitch(ji)   )
719               ENDIF
720               ! END MV MP 2016
721            END DO
722
723            ! Transfer ice energy to category jl2 by ridging
724            DO jk = 1, nlay_i
725               DO ij = 1, icells
726                  ji = indxi(ij) ; jj = indxj(ij)
727                  e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + erdg2(ij,jk) * fvol(ij) + eirft(ij,jk) * zswitch(ij)                 
728               END DO
729            END DO
730            !
731         END DO ! jl2
732         !
733      END DO ! jl1 (deforming categories)
734      !
735   END SUBROUTINE ice_rdgrft_ridgeshift
736
737
738   SUBROUTINE ice_rdgrft_icestrength( kstrngth )
739      !!----------------------------------------------------------------------
740      !!                ***  ROUTINE ice_rdgrft_icestrength ***
741      !!
742      !! ** Purpose :   computes ice strength used in dynamics routines of ice thickness
743      !!
744      !! ** Method  :   Compute the strength of the ice pack, defined as the energy (J m-2)
745      !!              dissipated per unit area removed from the ice pack under compression,
746      !!              and assumed proportional to the change in potential energy caused
747      !!              by ridging. Note that only Hibler's formulation is stable and that
748      !!              ice strength has to be smoothed
749      !!
750      !! ** Inputs / Ouputs : kstrngth (what kind of ice strength we are using)
751      !!----------------------------------------------------------------------
752      INTEGER, INTENT(in) ::   kstrngth    ! = 1 for Rothrock formulation, 0 for Hibler (1979)
753      !
754      INTEGER             ::   ji,jj, jl   ! dummy loop indices
755      INTEGER             ::   ismooth     ! smoothing the resistance to deformation
756      INTEGER             ::   itframe     ! number of time steps for the P smoothing
757      REAL(wp)            ::   zp, z1_3    ! local scalars
758      REAL(wp), DIMENSION(jpi,jpj) ::   zworka           ! temporary array used here
759      REAL(wp), DIMENSION(jpi,jpj) ::   zstrp1, zstrp2   ! strength at previous time steps
760      !!----------------------------------------------------------------------
761
762      !                              !--------------------------------------------------!
763      CALL ice_rdgrft_ridgeprep      ! Thickness distribution of ridging and ridged ice !
764      !                              !--------------------------------------------------!
765
766      !                              !--------------------------------------------------!
767      IF( kstrngth == 1 ) THEN       ! Ice strength => Rothrock (1975) method           !
768      !                              !--------------------------------------------------!
769         z1_3 = 1._wp / 3._wp
770         DO jl = 1, jpl
771            WHERE( athorn(:,:,jl) > 0._wp )
772               strength(:,:) =  -         athorn(:,:,jl) * ht_i(:,:,jl) * ht_i(:,:,jl)   &  ! PE loss from deforming ice
773                  &             + 2._wp * araft (:,:,jl) * ht_i(:,:,jl) * ht_i(:,:,jl)   &  ! PE gain from rafting ice
774                  &             +         aridge(:,:,jl) * krdg(:,:,jl) * z1_3 *   &        ! PE gain from ridging ice
775                  &                      ( hrmax(:,:,jl) * hrmax(:,:,jl) +         &
776                  &                        hrmin(:,:,jl) * hrmin(:,:,jl) +         &
777                  &                        hrmax(:,:,jl) * hrmin(:,:,jl) )
778            ELSEWHERE
779               strength(:,:) = 0._wp
780            END WHERE
781         END DO
782         strength(:,:) = rn_pe_rdg * zdrho * strength(:,:) / aksum(:,:) * tmask(:,:,1)
783                         ! where zdrho = (g/2)*(rhow-rhoi)*(rhoi/rhow) and rn_pe_rdg accounts for frictional dissipation
784         ismooth = 1
785      !                              !--------------------------------------------------!
786      ELSE                           ! Ice strength => Hibler (1979) method             !
787      !                              !--------------------------------------------------!
788         strength(:,:) = rn_pstar * vt_i(:,:) * EXP( - rn_crhg * ( 1._wp - at_i(:,:) )  ) * tmask(:,:,1)
789         !
790         ismooth = 1
791         !
792      ENDIF
793      !                              !--------------------------------------------------!
794      IF( ln_icestr_bvf ) THEN       ! Impact of brine volume                           !
795      !                              !--------------------------------------------------!
796      ! CAN BE REMOVED
797         DO jj = 1, jpj
798            DO ji = 1, jpi
799               strength(ji,jj) = strength(ji,jj) * exp(-5.88*SQRT(MAX(bvm_i(ji,jj),0.0)))
800            END DO
801         END DO
802      ENDIF
803      !                              !--------------------------------------------------!
804      SELECT CASE( ismooth )         ! Smoothing ice strength                           !
805      !                              !--------------------------------------------------!
806      !
807      CASE( 1 )               !--- Spatial smoothing
808         DO jj = 2, jpjm1
809            DO ji = 2, jpim1
810               IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp ) THEN
811                  zworka(ji,jj) = ( 4.0 * strength(ji,jj)              &
812                     &                  + strength(ji-1,jj) * tmask(ji-1,jj,1) + strength(ji+1,jj) * tmask(ji+1,jj,1) & 
813                     &                  + strength(ji,jj-1) * tmask(ji,jj-1,1) + strength(ji,jj+1) * tmask(ji,jj+1,1) &
814                     &            ) / ( 4.0 + tmask(ji-1,jj,1) + tmask(ji+1,jj,1) + tmask(ji,jj-1,1) + tmask(ji,jj+1,1) )
815               ELSE
816                  zworka(ji,jj) = 0._wp
817               ENDIF
818            END DO
819         END DO
820         
821         DO jj = 2, jpjm1
822            DO ji = 2, jpim1
823               strength(ji,jj) = zworka(ji,jj)
824            END DO
825         END DO
826         CALL lbc_lnk( strength, 'T', 1. )
827         !
828      CASE( 2 )               !--- Temporal smoothing
829         IF ( kt_ice == nit000 ) THEN
830            zstrp1(:,:) = 0._wp
831            zstrp2(:,:) = 0._wp
832         ENDIF
833         !
834         DO jj = 2, jpjm1
835            DO ji = 2, jpim1
836               IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp ) THEN
837                  itframe = 1 ! number of time steps for the running mean
838                  IF ( zstrp1(ji,jj) > 0._wp ) itframe = itframe + 1
839                  IF ( zstrp2(ji,jj) > 0._wp ) itframe = itframe + 1
840                  zp = ( strength(ji,jj) + zstrp1(ji,jj) + zstrp2(ji,jj) ) / itframe
841                  zstrp2  (ji,jj) = zstrp1  (ji,jj)
842                  zstrp1  (ji,jj) = strength(ji,jj)
843                  strength(ji,jj) = zp
844               ENDIF
845            END DO
846         END DO
847         CALL lbc_lnk( strength, 'T', 1. )
848         !
849      END SELECT
850      !
851   END SUBROUTINE ice_rdgrft_icestrength
852
853
854   SUBROUTINE ice_rdgrft_init
855      !!-------------------------------------------------------------------
856      !!                  ***  ROUTINE ice_rdgrft_init ***
857      !!
858      !! ** Purpose :   Physical constants and parameters linked
859      !!                to the mechanical ice redistribution
860      !!
861      !! ** Method  :   Read the namiceitdme namelist
862      !!                and check the parameters values
863      !!                called at the first timestep (nit000)
864      !!
865      !! ** input   :   Namelist namiceitdme
866      !!-------------------------------------------------------------------
867      INTEGER :: ios                 ! Local integer output status for namelist read
868      !!
869      NAMELIST/namiceitdme/ rn_cs     , nn_partfun, rn_gstar  , rn_astar   ,                & 
870         &                  ln_ridging, rn_hstar  , rn_por_rdg, rn_fsnowrdg, rn_fpondrdg,   & 
871         &                  ln_rafting, rn_hraft  , rn_craft  , rn_fsnowrft, rn_fpondrft
872      !!-------------------------------------------------------------------
873      !
874      REWIND( numnam_ice_ref )              ! Namelist namicetdme in reference namelist : Ice mechanical ice redistribution
875      READ  ( numnam_ice_ref, namiceitdme, IOSTAT = ios, ERR = 901)
876901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitdme in reference namelist', lwp )
877      !
878      REWIND( numnam_ice_cfg )              ! Namelist namiceitdme in configuration namelist : Ice mechanical ice redistribution
879      READ  ( numnam_ice_cfg, namiceitdme, IOSTAT = ios, ERR = 902 )
880902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitdme in configuration namelist', lwp )
881      IF(lwm) WRITE ( numoni, namiceitdme )
882      !
883      IF (lwp) THEN                          ! control print
884         WRITE(numout,*)
885         WRITE(numout,*) 'ice_rdgrft_init : ice parameters for mechanical ice redistribution '
886         WRITE(numout,*) '~~~~~~~~~~~~~~~'
887         WRITE(numout,*) '   Namelist namiceitdme'
888         WRITE(numout,*) '      Fraction of shear energy contributing to ridging        rn_cs       = ', rn_cs 
889         WRITE(numout,*) '      Switch for part. function (0) linear (1) exponential    nn_partfun  = ', nn_partfun
890         WRITE(numout,*) '      Fraction of total ice coverage contributing to ridging  rn_gstar    = ', rn_gstar
891         WRITE(numout,*) '      Equivalent to G* for an exponential part function       rn_astar    = ', rn_astar
892         WRITE(numout,*) '      Ridging of ice sheets or not                            ln_ridging  = ', ln_ridging
893         WRITE(numout,*) '      Quantity playing a role in max ridged ice thickness     rn_hstar    = ', rn_hstar
894         WRITE(numout,*) '      Initial porosity of ridges                              rn_por_rdg  = ', rn_por_rdg
895         WRITE(numout,*) '      Fraction of snow volume conserved during ridging        rn_fsnowrdg = ', rn_fsnowrdg 
896         WRITE(numout,*) '      Fraction of pond volume conserved during ridging        rn_fpondrdg = ', rn_fpondrdg 
897         WRITE(numout,*) '      Rafting of ice sheets or not                            ln_rafting  = ', ln_rafting
898         WRITE(numout,*) '      Parmeter thickness (threshold between ridge-raft)       rn_hraft    = ', rn_hraft
899         WRITE(numout,*) '      Rafting hyperbolic tangent coefficient                  rn_craft    = ', rn_craft 
900         WRITE(numout,*) '      Fraction of snow volume conserved during ridging        rn_fsnowrft = ', rn_fsnowrft 
901         WRITE(numout,*) '      Fraction of pond volume conserved during rafting        rn_fpondrft = ', rn_fpondrft 
902      ENDIF
903      !
904   END SUBROUTINE ice_rdgrft_init
905
906#else
907   !!----------------------------------------------------------------------
908   !!   Default option         Empty module          NO LIM-3 sea-ice model
909   !!----------------------------------------------------------------------
910#endif
911
912   !!======================================================================
913END MODULE icerdgrft
Note: See TracBrowser for help on using the repository browser.