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

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

source: branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 @ 6796

Last change on this file since 6796 was 6796, checked in by clem, 8 years ago

rewriting of ice rheology to conserve energy

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