New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
icerdgrft.F90 in branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

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

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

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