source: NEMO/releases/release-3.6/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90

Last change on this file was 7814, checked in by clem, 4 years ago

cosmetic changes for better phasing with trunk

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