source: branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 @ 7256

Last change on this file since 7256 was 7256, checked in by cbricaud, 4 years ago

phaze NEMO routines in CRS branch with nemo_v3_6_STABLE branch at rev 7213 (09-09-2016) (merge -r 5519:7213 )

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