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

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

changes in style - part5 - reaching the end

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