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.
limitd_me.F90 in branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 @ 8373

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

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