source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 3 years ago

All wrk_alloc removed

  • Property svn:keywords set to Id
File size: 45.6 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          ! 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
24   USE in_out_manager   ! I/O manager
25   USE iom              ! I/O manager
26   USE lib_fortran      ! glob_sum
27   USE timing           ! Timing
28   USE limcons          ! conservation tests
29   USE limctl           ! control prints
30
31   IMPLICIT NONE
32   PRIVATE
33
34   PUBLIC   lim_itd_me               ! called by ice_stp
35   PUBLIC   lim_itd_me_icestrength
36   PUBLIC   lim_itd_me_init
37   PUBLIC   lim_itd_me_alloc        ! called by sbc_lim_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$
61   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
62   !!----------------------------------------------------------------------
63CONTAINS
64
65   INTEGER FUNCTION lim_itd_me_alloc()
66      !!---------------------------------------------------------------------!
67      !!                ***  ROUTINE lim_itd_me_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=lim_itd_me_alloc )
74         !
75      IF( lim_itd_me_alloc /= 0 )   CALL ctl_warn( 'lim_itd_me_alloc: failed to allocate arrays' )
76      !
77   END FUNCTION lim_itd_me_alloc
78
79
80   SUBROUTINE lim_itd_me
81      !!---------------------------------------------------------------------!
82      !!                ***  ROUTINE lim_itd_me ***
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('limitd_me')
122
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      !
290      IF( nn_timing == 1 )  CALL timing_stop('limitd_me')
291   END SUBROUTINE lim_itd_me
292
293   SUBROUTINE lim_itd_me_ridgeprep
294      !!---------------------------------------------------------------------!
295      !!                ***  ROUTINE lim_itd_me_ridgeprep ***
296      !!
297      !! ** Purpose :   preparation for ridging and strength calculations
298      !!
299      !! ** Method  :   Compute the thickness distribution of the ice and open water
300      !!              participating in ridging and of the resulting ridges.
301      !!---------------------------------------------------------------------!
302      INTEGER ::   ji,jj, jl    ! dummy loop indices
303      REAL(wp) ::   Gstari, astari, hrmean, zdummy   ! local scalar
304      REAL(wp), DIMENSION(jpi,jpj,-1:jpl) ::   Gsum      ! Gsum(n) = sum of areas in categories 0 to n
305      !------------------------------------------------------------------------------!
306
307
308      Gstari     = 1.0/rn_gstar   
309      astari     = 1.0/rn_astar   
310      aksum(:,:)    = 0.0
311      athorn(:,:,:) = 0.0
312      aridge(:,:,:) = 0.0
313      araft (:,:,:) = 0.0
314
315      ! Zero out categories with very small areas
316      CALL lim_var_zapsmall
317
318      ! Ice thickness needed for rafting
319      DO jl = 1, jpl
320         DO jj = 1, jpj
321            DO ji = 1, jpi
322               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) )
323               ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
324           END DO
325         END DO
326      END DO
327
328      !------------------------------------------------------------------------------!
329      ! 1) Participation function
330      !------------------------------------------------------------------------------!
331
332      ! Compute total area of ice plus open water.
333      ! This is in general not equal to one because of divergence during transport
334      asum(:,:) = ato_i(:,:) + SUM( a_i, dim=3 )
335
336      ! Compute cumulative thickness distribution function
337      ! Compute the cumulative thickness distribution function Gsum,
338      ! where Gsum(n) is the fractional area in categories 0 to n.
339      ! initial value (in h = 0) equals open water area
340      Gsum(:,:,-1) = 0._wp
341      Gsum(:,:,0 ) = ato_i(:,:)
342      ! for each value of h, you have to add ice concentration then
343      DO jl = 1, jpl
344         Gsum(:,:,jl) = Gsum(:,:,jl-1) + a_i(:,:,jl)
345      END DO
346
347      ! Normalize the cumulative distribution to 1
348      DO jl = 0, jpl
349         Gsum(:,:,jl) = Gsum(:,:,jl) / asum(:,:)
350      END DO
351
352      ! 1.3 Compute participation function a(h) = b(h).g(h) (athorn)
353      !--------------------------------------------------------------------------------------------------
354      ! Compute the participation function athorn; this is analogous to
355      ! a(h) = b(h)g(h) as defined in Thorndike et al. (1975).
356      ! area lost from category n due to ridging/closing
357      ! athorn(n)   = total area lost due to ridging/closing
358      ! assume b(h) = (2/Gstar) * (1 - G(h)/Gstar).
359      !
360      ! The expressions for athorn are found by integrating b(h)g(h) between
361      ! the category boundaries.
362      ! athorn is always >= 0 and SUM(athorn(0:jpl))=1
363      !-----------------------------------------------------------------
364
365      IF( nn_partfun == 0 ) THEN       !--- Linear formulation (Thorndike et al., 1975)
366         DO jl = 0, jpl   
367            DO jj = 1, jpj 
368               DO ji = 1, jpi
369                  IF    ( Gsum(ji,jj,jl)   < rn_gstar ) THEN
370                     athorn(ji,jj,jl) = Gstari * ( Gsum(ji,jj,jl) - Gsum(ji,jj,jl-1) ) * &
371                        &                        ( 2._wp - ( Gsum(ji,jj,jl-1) + Gsum(ji,jj,jl) ) * Gstari )
372                  ELSEIF( Gsum(ji,jj,jl-1) < rn_gstar ) THEN
373                     athorn(ji,jj,jl) = Gstari * ( rn_gstar       - Gsum(ji,jj,jl-1) ) *  &
374                        &                        ( 2._wp - ( Gsum(ji,jj,jl-1) + rn_gstar       ) * Gstari )
375                  ELSE
376                     athorn(ji,jj,jl) = 0._wp
377                  ENDIF
378               END DO
379            END DO
380         END DO
381
382      ELSE                             !--- Exponential, more stable formulation (Lipscomb et al, 2007)
383         !                       
384         zdummy = 1._wp / ( 1._wp - EXP(-astari) )        ! precompute exponential terms using Gsum as a work array
385         DO jl = -1, jpl
386            Gsum(:,:,jl) = EXP( -Gsum(:,:,jl) * astari ) * zdummy
387         END DO
388         DO jl = 0, jpl
389             athorn(:,:,jl) = Gsum(:,:,jl-1) - Gsum(:,:,jl)
390         END DO
391         !
392      ENDIF
393
394      ! --- Ridging and rafting participation concentrations --- !
395      IF( ln_rafting .AND. ln_ridging ) THEN
396         !
397         DO jl = 1, jpl
398            DO jj = 1, jpj 
399               DO ji = 1, jpi
400                  zdummy           = TANH ( rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) )
401                  aridge(ji,jj,jl) = ( 1._wp + zdummy ) * 0.5_wp * athorn(ji,jj,jl)
402                  araft (ji,jj,jl) = athorn(ji,jj,jl) - aridge(ji,jj,jl)
403               END DO
404            END DO
405         END DO
406         !
407      ELSEIF( ln_ridging .AND. .NOT. ln_rafting ) THEN
408         !
409         DO jl = 1, jpl
410            aridge(:,:,jl) = athorn(:,:,jl)
411         END DO
412         !
413      ELSEIF( ln_rafting .AND. .NOT. ln_ridging ) THEN
414         !
415         DO jl = 1, jpl
416            araft(:,:,jl) = athorn(:,:,jl)
417         END DO
418         !
419      ENDIF
420
421      !-----------------------------------------------------------------
422      ! 2) Transfer function
423      !-----------------------------------------------------------------
424      ! Compute max and min ridged ice thickness for each ridging category.
425      ! Assume ridged ice is uniformly distributed between hrmin and hrmax.
426      !
427      ! This parameterization is a modified version of Hibler (1980).
428      ! The mean ridging thickness, hrmean, is proportional to hi^(0.5)
429      !  and for very thick ridging ice must be >= krdgmin*hi
430      !
431      ! The minimum ridging thickness, hrmin, is equal to 2*hi
432      !  (i.e., rafting) and for very thick ridging ice is
433      !  constrained by hrmin <= (hrmean + hi)/2.
434      !
435      ! The maximum ridging thickness, hrmax, is determined by
436      !  hrmean and hrmin.
437      !
438      ! These modifications have the effect of reducing the ice strength
439      ! (relative to the Hibler formulation) when very thick ice is
440      ! ridging.
441      !
442      ! aksum = net area removed/ total area removed
443      ! where total area removed = area of ice that ridges
444      !         net area removed = total area removed - area of new ridges
445      !-----------------------------------------------------------------
446
447      aksum(:,:) = athorn(:,:,0)
448      ! Transfer function
449      DO jl = 1, jpl !all categories have a specific transfer function
450         DO jj = 1, jpj
451            DO ji = 1, jpi
452               
453               IF( athorn(ji,jj,jl) > 0._wp ) THEN
454                  hrmean          = MAX( SQRT( rn_hstar * ht_i(ji,jj,jl) ), ht_i(ji,jj,jl) * krdgmin )
455                  hrmin(ji,jj,jl) = MIN( 2._wp * ht_i(ji,jj,jl), 0.5_wp * ( hrmean + ht_i(ji,jj,jl) ) )
456                  hrmax(ji,jj,jl) = 2._wp * hrmean - hrmin(ji,jj,jl)
457                  hraft(ji,jj,jl) = ht_i(ji,jj,jl) / kraft
458                  krdg(ji,jj,jl)  = ht_i(ji,jj,jl) / MAX( hrmean, epsi20 )
459
460                  ! Normalization factor : aksum, ensures mass conservation
461                  aksum(ji,jj) = aksum(ji,jj) + aridge(ji,jj,jl) * ( 1._wp - krdg(ji,jj,jl) )    &
462                     &                        + araft (ji,jj,jl) * ( 1._wp - kraft          )
463
464               ELSE
465                  hrmin(ji,jj,jl)  = 0._wp 
466                  hrmax(ji,jj,jl)  = 0._wp 
467                  hraft(ji,jj,jl)  = 0._wp 
468                  krdg (ji,jj,jl)  = 1._wp
469               ENDIF
470
471            END DO
472         END DO
473      END DO
474      !
475      !
476   END SUBROUTINE lim_itd_me_ridgeprep
477
478
479   SUBROUTINE lim_itd_me_ridgeshift( opning, closing_gross )
480      !!----------------------------------------------------------------------
481      !!                ***  ROUTINE lim_itd_me_icestrength ***
482      !!
483      !! ** Purpose :   shift ridging ice among thickness categories of ice thickness
484      !!
485      !! ** Method  :   Remove area, volume, and energy from each ridging category
486      !!              and add to thicker ice categories.
487      !!----------------------------------------------------------------------
488      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   opning         ! rate of opening due to divergence/shear
489      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   closing_gross  ! rate at which area removed, excluding area of new ridges
490      !
491      CHARACTER (len=80) ::   fieldid   ! field identifier
492      !
493      INTEGER ::   ji, jj, jl, jl1, jl2, jk   ! dummy loop indices
494      INTEGER ::   ij                ! horizontal index, combines i and j loops
495      INTEGER ::   icells            ! number of cells with a_i > puny
496      REAL(wp) ::   hL, hR, farea    ! left and right limits of integration
497
498      INTEGER , DIMENSION(jpij) ::   indxi, indxj   ! compressed indices
499      REAL(wp), DIMENSION(jpij) ::   zswitch, fvol   ! new ridge volume going to n2
500
501      REAL(wp), DIMENSION(jpij) ::   afrac            ! fraction of category area ridged
502      REAL(wp), DIMENSION(jpij) ::   ardg1 , ardg2    ! area of ice ridged & new ridges
503      REAL(wp), DIMENSION(jpij) ::   vsrdg , esrdg    ! snow volume & energy of ridging ice
504      REAL(wp), DIMENSION(jpij) ::   dhr   , dhr2     ! hrmax - hrmin  &  hrmax^2 - hrmin^2
505
506      REAL(wp), DIMENSION(jpij) ::   vrdg1   ! volume of ice ridged
507      REAL(wp), DIMENSION(jpij) ::   vrdg2   ! volume of new ridges
508      REAL(wp), DIMENSION(jpij) ::   vsw     ! volume of seawater trapped into ridges
509      REAL(wp), DIMENSION(jpij) ::   srdg1   ! sal*volume of ice ridged
510      REAL(wp), DIMENSION(jpij) ::   srdg2   ! sal*volume of new ridges
511      REAL(wp), DIMENSION(jpij) ::   smsw    ! sal*volume of water trapped into ridges
512      REAL(wp), DIMENSION(jpij) ::   oirdg1, oirdg2   ! ice age of ice ridged
513
514      REAL(wp), DIMENSION(jpij) ::   afrft            ! fraction of category area rafted
515      REAL(wp), DIMENSION(jpij) ::   arft1 , arft2    ! area of ice rafted and new rafted zone
516      REAL(wp), DIMENSION(jpij) ::   virft , vsrft    ! ice & snow volume of rafting ice
517      REAL(wp), DIMENSION(jpij) ::   esrft , smrft    ! snow energy & salinity of rafting ice
518      REAL(wp), DIMENSION(jpij) ::   oirft1, oirft2   ! ice age of ice rafted
519
520      REAL(wp), DIMENSION(jpij,nlay_i) ::   eirft      ! ice energy of rafting ice
521      REAL(wp), DIMENSION(jpij,nlay_i) ::   erdg1      ! enth*volume of ice ridged
522      REAL(wp), DIMENSION(jpij,nlay_i) ::   erdg2      ! enth*volume of new ridges
523      REAL(wp), DIMENSION(jpij,nlay_i) ::   ersw       ! enth of water trapped into ridges
524      !!----------------------------------------------------------------------
525
526
527      !-------------------------------------------------------------------------------
528      ! 1) Compute change in open water area due to closing and opening.
529      !-------------------------------------------------------------------------------
530      DO jj = 1, jpj
531         DO ji = 1, jpi
532            ato_i(ji,jj) = MAX( 0._wp, ato_i(ji,jj) +  &
533               &                     ( opning(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) ) * rdt_ice )
534         END DO
535      END DO
536
537      !-----------------------------------------------------------------
538      ! 3) Pump everything from ice which is being ridged / rafted
539      !-----------------------------------------------------------------
540      ! Compute the area, volume, and energy of ice ridging in each
541      ! category, along with the area of the resulting ridge.
542
543      DO jl1 = 1, jpl !jl1 describes the ridging category
544
545         !------------------------------------------------
546         ! 3.1) Identify grid cells with nonzero ridging
547         !------------------------------------------------
548         icells = 0
549         DO jj = 1, jpj
550            DO ji = 1, jpi
551               IF( athorn(ji,jj,jl1) > 0._wp .AND. closing_gross(ji,jj) > 0._wp ) THEN
552                  icells = icells + 1
553                  indxi(icells) = ji
554                  indxj(icells) = jj
555               ENDIF
556            END DO
557         END DO
558
559         DO ij = 1, icells
560            ji = indxi(ij) ; jj = indxj(ij)
561
562            !--------------------------------------------------------------------
563            ! 3.2) Compute area of ridging ice (ardg1) and of new ridge (ardg2)
564            !--------------------------------------------------------------------
565            ardg1(ij) = aridge(ji,jj,jl1) * closing_gross(ji,jj) * rdt_ice
566            arft1(ij) = araft (ji,jj,jl1) * closing_gross(ji,jj) * rdt_ice
567
568            !---------------------------------------------------------------
569            ! 3.3) Compute ridging /rafting fractions, make sure afrac <=1
570            !---------------------------------------------------------------
571            afrac(ij) = ardg1(ij) / a_i(ji,jj,jl1) !ridging
572            afrft(ij) = arft1(ij) / a_i(ji,jj,jl1) !rafting
573            ardg2(ij) = ardg1(ij) * krdg(ji,jj,jl1)
574            arft2(ij) = arft1(ij) * kraft
575
576            !--------------------------------------------------------------------------
577            ! 3.4) Subtract area, volume, and energy from ridging
578            !     / rafting category n1.
579            !--------------------------------------------------------------------------
580            vrdg1(ij) = v_i(ji,jj,jl1) * afrac(ij)
581            vrdg2(ij) = vrdg1(ij) * ( 1. + rn_por_rdg )
582            vsw  (ij) = vrdg1(ij) * rn_por_rdg
583
584            vsrdg (ij) = v_s  (ji,jj,  jl1) * afrac(ij)
585            esrdg (ij) = e_s  (ji,jj,1,jl1) * afrac(ij)
586            srdg1 (ij) = smv_i(ji,jj,  jl1) * afrac(ij)
587            oirdg1(ij) = oa_i (ji,jj,  jl1) * afrac(ij)
588            oirdg2(ij) = oa_i (ji,jj,  jl1) * afrac(ij) * krdg(ji,jj,jl1) 
589
590            ! rafting volumes, heat contents ...
591            virft (ij) = v_i  (ji,jj,  jl1) * afrft(ij)
592            vsrft (ij) = v_s  (ji,jj,  jl1) * afrft(ij)
593            esrft (ij) = e_s  (ji,jj,1,jl1) * afrft(ij)
594            smrft (ij) = smv_i(ji,jj,  jl1) * afrft(ij) 
595            oirft1(ij) = oa_i (ji,jj,  jl1) * afrft(ij) 
596            oirft2(ij) = oa_i (ji,jj,  jl1) * afrft(ij) * kraft 
597
598            !-----------------------------------------------------------------
599            ! 3.5) Compute properties of new ridges
600            !-----------------------------------------------------------------
601            smsw(ij)  = vsw(ij) * sss_m(ji,jj)                   ! salt content of seawater frozen in voids
602            srdg2(ij) = srdg1(ij) + smsw(ij)                     ! salt content of new ridge
603           
604            sfx_dyn(ji,jj) = sfx_dyn(ji,jj) - smsw(ij) * rhoic * r1_rdtice
605            wfx_dyn(ji,jj) = wfx_dyn(ji,jj) - vsw (ij) * rhoic * r1_rdtice   ! increase in ice volume due to seawater frozen in voids
606
607            ! virtual salt flux to keep salinity constant
608            IF( nn_icesal == 1 .OR. nn_icesal == 3 )  THEN
609               srdg2(ij)      = srdg2(ij) - vsw(ij) * ( sss_m(ji,jj) - sm_i(ji,jj,jl1) )           ! ridge salinity = sm_i
610               sfx_bri(ji,jj) = sfx_bri(ji,jj) + sss_m(ji,jj)    * vsw(ij) * rhoic * r1_rdtice  &  ! put back sss_m into the ocean
611                  &                            - sm_i(ji,jj,jl1) * vsw(ij) * rhoic * r1_rdtice     ! and get  sm_i  from the ocean
612            ENDIF
613               
614            !------------------------------------------           
615            ! 3.7 Put the snow somewhere in the ocean
616            !------------------------------------------           
617            !  Place part of the snow lost by ridging into the ocean.
618            !  Note that esrdg > 0; the ocean must cool to melt snow.
619            !  If the ocean temp = Tf already, new ice must grow.
620            !  During the next time step, thermo_rates will determine whether
621            !  the ocean cools or new ice grows.
622            wfx_snw(ji,jj) = wfx_snw(ji,jj) + ( rhosn * vsrdg(ij) * ( 1._wp - rn_fsnowrdg )   & 
623               &                              + rhosn * vsrft(ij) * ( 1._wp - rn_fsnowrft ) ) * r1_rdtice  ! fresh water source for ocean
624
625            hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ( - esrdg(ij) * ( 1._wp - rn_fsnowrdg )         & 
626               &                                - esrft(ij) * ( 1._wp - rn_fsnowrft ) ) * r1_rdtice        ! heat sink for ocean (<0, W.m-2)
627
628            !-----------------------------------------------------------------
629            ! 3.8 Compute quantities used to apportion ice among categories
630            ! in the n2 loop below
631            !-----------------------------------------------------------------
632            dhr (ij) = 1._wp / ( hrmax(ji,jj,jl1)                    - hrmin(ji,jj,jl1)                    )
633            dhr2(ij) = 1._wp / ( hrmax(ji,jj,jl1) * hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) * hrmin(ji,jj,jl1) )
634
635
636            ! update jl1 (removing ridged/rafted area)
637            a_i  (ji,jj,  jl1) = a_i  (ji,jj,  jl1) - ardg1 (ij) - arft1 (ij)
638            v_i  (ji,jj,  jl1) = v_i  (ji,jj,  jl1) - vrdg1 (ij) - virft (ij)
639            v_s  (ji,jj,  jl1) = v_s  (ji,jj,  jl1) - vsrdg (ij) - vsrft (ij)
640            e_s  (ji,jj,1,jl1) = e_s  (ji,jj,1,jl1) - esrdg (ij) - esrft (ij)
641            smv_i(ji,jj,  jl1) = smv_i(ji,jj,  jl1) - srdg1 (ij) - smrft (ij)
642            oa_i (ji,jj,  jl1) = oa_i (ji,jj,  jl1) - oirdg1(ij) - oirft1(ij)
643
644         END DO
645
646         !--------------------------------------------------------------------
647         ! 3.9 Compute ridging ice enthalpy, remove it from ridging ice and
648         !      compute ridged ice enthalpy
649         !--------------------------------------------------------------------
650         DO jk = 1, nlay_i
651            DO ij = 1, icells
652               ji = indxi(ij) ; jj = indxj(ij)
653               ! heat content of ridged ice
654               erdg1(ij,jk) = e_i(ji,jj,jk,jl1) * afrac(ij) 
655               eirft(ij,jk) = e_i(ji,jj,jk,jl1) * afrft(ij)               
656               
657               ! enthalpy of the trapped seawater (J/m2, >0)
658               ! clem: if sst>0, then ersw <0 (is that possible?)
659               ersw(ij,jk)  = - rhoic * vsw(ij) * rcp * sst_m(ji,jj) * r1_nlay_i
660
661               ! heat flux to the ocean
662               hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ersw(ij,jk) * r1_rdtice  ! > 0 [W.m-2] ocean->ice flux
663
664               ! it is added to sea ice because the sign convention is the opposite of the sign convention for the ocean
665               erdg2(ij,jk) = erdg1(ij,jk) + ersw(ij,jk)
666
667               ! update jl1
668               e_i  (ji,jj,jk,jl1) = e_i(ji,jj,jk,jl1) - erdg1(ij,jk) - eirft(ij,jk)
669
670            END DO
671         END DO
672
673         !-------------------------------------------------------------------------------
674         ! 4) Add area, volume, and energy of new ridge to each category jl2
675         !-------------------------------------------------------------------------------
676         DO jl2  = 1, jpl 
677            ! over categories to which ridged/rafted ice is transferred
678            DO ij = 1, icells
679               ji = indxi(ij) ; jj = indxj(ij)
680
681               ! Compute the fraction of ridged ice area and volume going to thickness category jl2.
682               IF( hrmin(ji,jj,jl1) <= hi_max(jl2) .AND. hrmax(ji,jj,jl1) > hi_max(jl2-1) ) THEN
683                  hL = MAX( hrmin(ji,jj,jl1), hi_max(jl2-1) )
684                  hR = MIN( hrmax(ji,jj,jl1), hi_max(jl2)   )
685                  farea    = ( hR      - hL      ) * dhr(ij) 
686                  fvol(ij) = ( hR * hR - hL * hL ) * dhr2(ij)
687               ELSE
688                  farea    = 0._wp 
689                  fvol(ij) = 0._wp                 
690               ENDIF
691
692               ! Compute the fraction of rafted ice area and volume going to thickness category jl2
693               IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. hraft(ji,jj,jl1) >  hi_max(jl2-1) ) THEN
694                  zswitch(ij) = 1._wp
695               ELSE
696                  zswitch(ij) = 0._wp                 
697               ENDIF
698
699               a_i  (ji,jj  ,jl2) = a_i  (ji,jj  ,jl2) + ( ardg2 (ij) * farea    + arft2 (ij) * zswitch(ij) )
700               oa_i (ji,jj  ,jl2) = oa_i (ji,jj  ,jl2) + ( oirdg2(ij) * farea    + oirft2(ij) * zswitch(ij) )
701               v_i  (ji,jj  ,jl2) = v_i  (ji,jj  ,jl2) + ( vrdg2 (ij) * fvol(ij) + virft (ij) * zswitch(ij) )
702               smv_i(ji,jj  ,jl2) = smv_i(ji,jj  ,jl2) + ( srdg2 (ij) * fvol(ij) + smrft (ij) * zswitch(ij) )
703               v_s  (ji,jj  ,jl2) = v_s  (ji,jj  ,jl2) + ( vsrdg (ij) * rn_fsnowrdg * fvol(ij)  +  &
704                  &                                        vsrft (ij) * rn_fsnowrft * zswitch(ij) )
705               e_s  (ji,jj,1,jl2) = e_s  (ji,jj,1,jl2) + ( esrdg (ij) * rn_fsnowrdg * fvol(ij)  +  &
706                  &                                        esrft (ij) * rn_fsnowrft * zswitch(ij) )
707
708            END DO
709
710            ! Transfer ice energy to category jl2 by ridging
711            DO jk = 1, nlay_i
712               DO ij = 1, icells
713                  ji = indxi(ij) ; jj = indxj(ij)
714                  e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + erdg2(ij,jk) * fvol(ij) + eirft(ij,jk) * zswitch(ij)                 
715               END DO
716            END DO
717            !
718         END DO ! jl2
719         
720      END DO ! jl1 (deforming categories)
721
722      !
723      !
724   END SUBROUTINE lim_itd_me_ridgeshift
725
726   SUBROUTINE lim_itd_me_icestrength( kstrngth )
727      !!----------------------------------------------------------------------
728      !!                ***  ROUTINE lim_itd_me_icestrength ***
729      !!
730      !! ** Purpose :   computes ice strength used in dynamics routines of ice thickness
731      !!
732      !! ** Method  :   Compute the strength of the ice pack, defined as the energy (J m-2)
733      !!              dissipated per unit area removed from the ice pack under compression,
734      !!              and assumed proportional to the change in potential energy caused
735      !!              by ridging. Note that only Hibler's formulation is stable and that
736      !!              ice strength has to be smoothed
737      !!
738      !! ** Inputs / Ouputs : kstrngth (what kind of ice strength we are using)
739      !!----------------------------------------------------------------------
740      INTEGER, INTENT(in) ::   kstrngth    ! = 1 for Rothrock formulation, 0 for Hibler (1979)
741      INTEGER             ::   ji,jj, jl   ! dummy loop indices
742      INTEGER             ::   ksmooth     ! smoothing the resistance to deformation
743      INTEGER             ::   numts_rm    ! number of time steps for the P smoothing
744      REAL(wp)            ::   zp, z1_3    ! local scalars
745      REAL(wp), DIMENSION(jpi,jpj) ::   zworka           ! temporary array used here
746      REAL(wp), DIMENSION(jpi,jpj) ::   zstrp1, zstrp2   ! strength at previous time steps
747      !!----------------------------------------------------------------------
748
749
750      !------------------------------------------------------------------------------!
751      ! 1) Initialize
752      !------------------------------------------------------------------------------!
753      strength(:,:) = 0._wp
754
755      !------------------------------------------------------------------------------!
756      ! 2) Compute thickness distribution of ridging and ridged ice
757      !------------------------------------------------------------------------------!
758      CALL lim_itd_me_ridgeprep
759
760      !------------------------------------------------------------------------------!
761      ! 3) Rothrock(1975)'s method
762      !------------------------------------------------------------------------------!
763      IF( kstrngth == 1 ) THEN
764         z1_3 = 1._wp / 3._wp
765         DO jl = 1, jpl
766            DO jj= 1, jpj
767               DO ji = 1, jpi
768                  !
769                  IF( athorn(ji,jj,jl) > 0._wp ) THEN
770                     !----------------------------
771                     ! PE loss from deforming ice
772                     !----------------------------
773                     strength(ji,jj) = strength(ji,jj) - athorn(ji,jj,jl) * ht_i(ji,jj,jl) * ht_i(ji,jj,jl)
774
775                     !--------------------------
776                     ! PE gain from rafting ice
777                     !--------------------------
778                     strength(ji,jj) = strength(ji,jj) + 2._wp * araft(ji,jj,jl) * ht_i(ji,jj,jl) * ht_i(ji,jj,jl)
779
780                     !----------------------------
781                     ! PE gain from ridging ice
782                     !----------------------------
783                     strength(ji,jj) = strength(ji,jj) + aridge(ji,jj,jl) * krdg(ji,jj,jl) * z1_3 *  &
784                        &                              ( hrmax(ji,jj,jl) * hrmax(ji,jj,jl) +         &
785                        &                                hrmin(ji,jj,jl) * hrmin(ji,jj,jl) +         &
786                        &                                hrmax(ji,jj,jl) * hrmin(ji,jj,jl) ) 
787                        !!(a**3-b**3)/(a-b) = a*a+ab+b*b                     
788                  ENDIF
789                  !
790               END DO
791            END DO
792         END DO
793   
794         strength(:,:) = rn_pe_rdg * Cp * strength(:,:) / aksum(:,:) * tmask(:,:,1)
795                         ! where Cp = (g/2)*(rhow-rhoi)*(rhoi/rhow) and rn_pe_rdg accounts for frictional dissipation
796         ksmooth = 1
797
798      !------------------------------------------------------------------------------!
799      ! 4) Hibler (1979)' method
800      !------------------------------------------------------------------------------!
801      ELSE                      ! kstrngth ne 1:  Hibler (1979) form
802         !
803         strength(:,:) = rn_pstar * vt_i(:,:) * EXP( - rn_crhg * ( 1._wp - at_i(:,:) )  ) * tmask(:,:,1)
804         !
805         ksmooth = 1
806         !
807      ENDIF                     ! kstrngth
808      !
809      !------------------------------------------------------------------------------!
810      ! 5) Impact of brine volume
811      !------------------------------------------------------------------------------!
812      ! CAN BE REMOVED
813      IF( ln_icestr_bvf ) THEN
814         DO jj = 1, jpj
815            DO ji = 1, jpi
816               strength(ji,jj) = strength(ji,jj) * exp(-5.88*SQRT(MAX(bvm_i(ji,jj),0.0)))
817            END DO
818         END DO
819      ENDIF
820      !
821      !------------------------------------------------------------------------------!
822      ! 6) Smoothing ice strength
823      !------------------------------------------------------------------------------!
824      !
825      !-------------------
826      ! Spatial smoothing
827      !-------------------
828      IF ( ksmooth == 1 ) THEN
829
830         DO jj = 2, jpjm1
831            DO ji = 2, jpim1
832               IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp ) THEN
833                  zworka(ji,jj) = ( 4.0 * strength(ji,jj)              &
834                     &                  + strength(ji-1,jj) * tmask(ji-1,jj,1) + strength(ji+1,jj) * tmask(ji+1,jj,1) & 
835                     &                  + strength(ji,jj-1) * tmask(ji,jj-1,1) + strength(ji,jj+1) * tmask(ji,jj+1,1) &
836                     &            ) / ( 4.0 + tmask(ji-1,jj,1) + tmask(ji+1,jj,1) + tmask(ji,jj-1,1) + tmask(ji,jj+1,1) )
837               ELSE
838                  zworka(ji,jj) = 0._wp
839               ENDIF
840            END DO
841         END DO
842
843         DO jj = 2, jpjm1
844            DO ji = 2, jpim1
845               strength(ji,jj) = zworka(ji,jj)
846            END DO
847         END DO
848         CALL lbc_lnk( strength, 'T', 1. )
849
850      ENDIF
851
852      !--------------------
853      ! Temporal smoothing
854      !--------------------
855      IF ( ksmooth == 2 ) THEN
856
857         IF ( numit == nit000 + nn_fsbc - 1 ) THEN
858            zstrp1(:,:) = 0._wp
859            zstrp2(:,:) = 0._wp
860         ENDIF
861
862         DO jj = 2, jpjm1
863            DO ji = 2, jpim1
864               IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp ) THEN
865                  numts_rm = 1 ! number of time steps for the running mean
866                  IF ( zstrp1(ji,jj) > 0._wp ) numts_rm = numts_rm + 1
867                  IF ( zstrp2(ji,jj) > 0._wp ) numts_rm = numts_rm + 1
868                  zp = ( strength(ji,jj) + zstrp1(ji,jj) + zstrp2(ji,jj) ) / numts_rm
869                  zstrp2(ji,jj) = zstrp1(ji,jj)
870                  zstrp1(ji,jj) = strength(ji,jj)
871                  strength(ji,jj) = zp
872               ENDIF
873            END DO
874         END DO
875
876         CALL lbc_lnk( strength, 'T', 1. )      ! Boundary conditions
877
878      ENDIF ! ksmooth
879
880      !
881   END SUBROUTINE lim_itd_me_icestrength
882
883   SUBROUTINE lim_itd_me_init
884      !!-------------------------------------------------------------------
885      !!                   ***  ROUTINE lim_itd_me_init ***
886      !!
887      !! ** Purpose :   Physical constants and parameters linked
888      !!                to the mechanical ice redistribution
889      !!
890      !! ** Method  :   Read the namiceitdme namelist
891      !!                and check the parameters values
892      !!                called at the first timestep (nit000)
893      !!
894      !! ** input   :   Namelist namiceitdme
895      !!-------------------------------------------------------------------
896      INTEGER :: ios                 ! Local integer output status for namelist read
897      NAMELIST/namiceitdme/ rn_cs, nn_partfun, rn_gstar, rn_astar,             & 
898        &                   ln_ridging, rn_hstar, rn_por_rdg, rn_fsnowrdg, ln_rafting, rn_hraft, rn_craft, rn_fsnowrft
899      !!-------------------------------------------------------------------
900      !
901      REWIND( numnam_ice_ref )              ! Namelist namicetdme in reference namelist : Ice mechanical ice redistribution
902      READ  ( numnam_ice_ref, namiceitdme, IOSTAT = ios, ERR = 901)
903901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitdme in reference namelist', lwp )
904
905      REWIND( numnam_ice_cfg )              ! Namelist namiceitdme in configuration namelist : Ice mechanical ice redistribution
906      READ  ( numnam_ice_cfg, namiceitdme, IOSTAT = ios, ERR = 902 )
907902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitdme in configuration namelist', lwp )
908      IF(lwm) WRITE ( numoni, namiceitdme )
909      !
910      IF (lwp) THEN                          ! control print
911         WRITE(numout,*)
912         WRITE(numout,*)'lim_itd_me_init : ice parameters for mechanical ice redistribution '
913         WRITE(numout,*)'~~~~~~~~~~~~~~~'
914         WRITE(numout,*)'   Fraction of shear energy contributing to ridging        rn_cs       = ', rn_cs 
915         WRITE(numout,*)'   Switch for part. function (0) linear (1) exponential    nn_partfun  = ', nn_partfun
916         WRITE(numout,*)'   Fraction of total ice coverage contributing to ridging  rn_gstar    = ', rn_gstar
917         WRITE(numout,*)'   Equivalent to G* for an exponential part function       rn_astar    = ', rn_astar
918         WRITE(numout,*)'   Ridging of ice sheets or not                            ln_ridging  = ', ln_ridging
919         WRITE(numout,*)'   Quantity playing a role in max ridged ice thickness     rn_hstar    = ', rn_hstar
920         WRITE(numout,*)'   Initial porosity of ridges                              rn_por_rdg  = ', rn_por_rdg
921         WRITE(numout,*)'   Fraction of snow volume conserved during ridging        rn_fsnowrdg = ', rn_fsnowrdg 
922         WRITE(numout,*)'   Rafting of ice sheets or not                            ln_rafting  = ', ln_rafting
923         WRITE(numout,*)'   Parmeter thickness (threshold between ridge-raft)       rn_hraft    = ', rn_hraft
924         WRITE(numout,*)'   Rafting hyperbolic tangent coefficient                  rn_craft    = ', rn_craft 
925         WRITE(numout,*)'   Fraction of snow volume conserved during ridging        rn_fsnowrft = ', rn_fsnowrft 
926      ENDIF
927      !
928   END SUBROUTINE lim_itd_me_init
929
930#else
931   !!----------------------------------------------------------------------
932   !!   Default option         Empty module          NO LIM-3 sea-ice model
933   !!----------------------------------------------------------------------
934CONTAINS
935   SUBROUTINE lim_itd_me           ! Empty routines
936   END SUBROUTINE lim_itd_me
937   SUBROUTINE lim_itd_me_icestrength
938   END SUBROUTINE lim_itd_me_icestrength
939   SUBROUTINE lim_itd_me_init
940   END SUBROUTINE lim_itd_me_init
941#endif
942   !!======================================================================
943END MODULE limitd_me
Note: See TracBrowser for help on using the repository browser.