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 @ 6584

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

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