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

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

changes in style - part6 - one more round

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