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/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 @ 6498

Last change on this file since 6498 was 6498, checked in by timgraham, 8 years ago

Merge head of nemo_v3_6_STABLE into package branch

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