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

source: branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 @ 8150

Last change on this file since 8150 was 8150, checked in by vancop, 7 years ago

SIMIP outputs, phase 2
  • Property svn:keywords set to Id
File size: 49.2 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      REAL(wp) ::   zwfx_snw         ! snow mass flux increment
538
539      INTEGER , POINTER, DIMENSION(:) ::   indxi, indxj   ! compressed indices
540      REAL(wp), POINTER, DIMENSION(:) ::   zswitch, fvol   ! new ridge volume going to n2
541
542      REAL(wp), POINTER, DIMENSION(:) ::   afrac            ! fraction of category area ridged
543      REAL(wp), POINTER, DIMENSION(:) ::   ardg1 , ardg2    ! area of ice ridged & new ridges
544      REAL(wp), POINTER, DIMENSION(:) ::   vsrdg , esrdg    ! snow volume & energy of ridging ice
545      REAL(wp), POINTER, DIMENSION(:) ::   dhr   , dhr2     ! hrmax - hrmin  &  hrmax^2 - hrmin^2
546
547      REAL(wp), POINTER, DIMENSION(:) ::   vrdg1   ! volume of ice ridged
548      REAL(wp), POINTER, DIMENSION(:) ::   vrdg2   ! volume of new ridges
549      REAL(wp), POINTER, DIMENSION(:) ::   vsw     ! volume of seawater trapped into ridges
550      REAL(wp), POINTER, DIMENSION(:) ::   srdg1   ! sal*volume of ice ridged
551      REAL(wp), POINTER, DIMENSION(:) ::   srdg2   ! sal*volume of new ridges
552      REAL(wp), POINTER, DIMENSION(:) ::   smsw    ! sal*volume of water trapped into ridges
553      REAL(wp), POINTER, DIMENSION(:) ::   oirdg1, oirdg2   ! ice age of ice ridged
554
555      REAL(wp), POINTER, DIMENSION(:) ::   afrft            ! fraction of category area rafted
556      REAL(wp), POINTER, DIMENSION(:) ::   arft1 , arft2    ! area of ice rafted and new rafted zone
557      REAL(wp), POINTER, DIMENSION(:) ::   virft , vsrft    ! ice & snow volume of rafting ice
558      REAL(wp), POINTER, DIMENSION(:) ::   esrft , smrft    ! snow energy & salinity of rafting ice
559      REAL(wp), POINTER, DIMENSION(:) ::   oirft1, oirft2   ! ice age of ice rafted
560
561      REAL(wp), POINTER, DIMENSION(:,:) ::   eirft      ! ice energy of rafting ice
562      REAL(wp), POINTER, DIMENSION(:,:) ::   erdg1      ! enth*volume of ice ridged
563      REAL(wp), POINTER, DIMENSION(:,:) ::   erdg2      ! enth*volume of new ridges
564      REAL(wp), POINTER, DIMENSION(:,:) ::   ersw       ! enth of water trapped into ridges
565      !!----------------------------------------------------------------------
566
567      CALL wrk_alloc( jpij,        indxi, indxj )
568      CALL wrk_alloc( jpij,        zswitch, fvol )
569      CALL wrk_alloc( jpij,        afrac, ardg1, ardg2, vsrdg, esrdg, dhr, dhr2 )
570      CALL wrk_alloc( jpij,        vrdg1, vrdg2, vsw  , srdg1, srdg2, smsw, oirdg1, oirdg2 )
571      CALL wrk_alloc( jpij,        afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 )
572      CALL wrk_alloc( jpij,nlay_i, eirft, erdg1, erdg2, ersw )
573
574      !-------------------------------------------------------------------------------
575      ! 1) Compute change in open water area due to closing and opening.
576      !-------------------------------------------------------------------------------
577      DO jj = 1, jpj
578         DO ji = 1, jpi
579            ato_i(ji,jj) = MAX( 0._wp, ato_i(ji,jj) +  &
580               &                     ( opning(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) ) * rdt_ice )
581         END DO
582      END DO
583
584      !-----------------------------------------------------------------
585      ! 3) Pump everything from ice which is being ridged / rafted
586      !-----------------------------------------------------------------
587      ! Compute the area, volume, and energy of ice ridging in each
588      ! category, along with the area of the resulting ridge.
589
590      DO jl1 = 1, jpl !jl1 describes the ridging category
591
592         !------------------------------------------------
593         ! 3.1) Identify grid cells with nonzero ridging
594         !------------------------------------------------
595         icells = 0
596         DO jj = 1, jpj
597            DO ji = 1, jpi
598               IF( athorn(ji,jj,jl1) > 0._wp .AND. closing_gross(ji,jj) > 0._wp ) THEN
599                  icells = icells + 1
600                  indxi(icells) = ji
601                  indxj(icells) = jj
602               ENDIF
603            END DO
604         END DO
605
606         DO ij = 1, icells
607            ji = indxi(ij) ; jj = indxj(ij)
608
609            !--------------------------------------------------------------------
610            ! 3.2) Compute area of ridging ice (ardg1) and of new ridge (ardg2)
611            !--------------------------------------------------------------------
612            ardg1(ij) = aridge(ji,jj,jl1) * closing_gross(ji,jj) * rdt_ice
613            arft1(ij) = araft (ji,jj,jl1) * closing_gross(ji,jj) * rdt_ice
614
615            !---------------------------------------------------------------
616            ! 3.3) Compute ridging /rafting fractions, make sure afrac <=1
617            !---------------------------------------------------------------
618            afrac(ij) = ardg1(ij) / a_i(ji,jj,jl1) !ridging
619            afrft(ij) = arft1(ij) / a_i(ji,jj,jl1) !rafting
620            ardg2(ij) = ardg1(ij) * krdg(ji,jj,jl1)
621            arft2(ij) = arft1(ij) * kraft
622
623            !--------------------------------------------------------------------------
624            ! 3.4) Subtract area, volume, and energy from ridging
625            !     / rafting category n1.
626            !--------------------------------------------------------------------------
627            vrdg1(ij) = v_i(ji,jj,jl1) * afrac(ij)
628            vrdg2(ij) = vrdg1(ij) * ( 1. + rn_por_rdg )
629            vsw  (ij) = vrdg1(ij) * rn_por_rdg
630
631            vsrdg (ij) = v_s  (ji,jj,  jl1) * afrac(ij)
632            esrdg (ij) = e_s  (ji,jj,1,jl1) * afrac(ij)
633            srdg1 (ij) = smv_i(ji,jj,  jl1) * afrac(ij)
634            oirdg1(ij) = oa_i (ji,jj,  jl1) * afrac(ij)
635            oirdg2(ij) = oa_i (ji,jj,  jl1) * afrac(ij) * krdg(ji,jj,jl1) 
636
637            ! rafting volumes, heat contents ...
638            virft (ij) = v_i  (ji,jj,  jl1) * afrft(ij)
639            vsrft (ij) = v_s  (ji,jj,  jl1) * afrft(ij)
640            esrft (ij) = e_s  (ji,jj,1,jl1) * afrft(ij)
641            smrft (ij) = smv_i(ji,jj,  jl1) * afrft(ij) 
642            oirft1(ij) = oa_i (ji,jj,  jl1) * afrft(ij) 
643            oirft2(ij) = oa_i (ji,jj,  jl1) * afrft(ij) * kraft 
644
645            !-----------------------------------------------------------------
646            ! 3.5) Compute properties of new ridges
647            !-----------------------------------------------------------------
648            smsw(ij)  = vsw(ij) * sss_m(ji,jj)                   ! salt content of seawater frozen in voids
649            srdg2(ij) = srdg1(ij) + smsw(ij)                     ! salt content of new ridge
650           
651            sfx_dyn(ji,jj)     = sfx_dyn(ji,jj)     - smsw(ij) * rhoic * r1_rdtice
652            wfx_dyn(ji,jj)     = wfx_dyn(ji,jj)     - vsw (ij) * rhoic * r1_rdtice   ! increase in ice volume due to seawater frozen in voids
653           
654             ! virtual salt flux to keep salinity constant
655            IF( nn_icesal == 1 .OR. nn_icesal == 3 )  THEN
656               srdg2(ij)      = srdg2(ij) - vsw(ij) * ( sss_m(ji,jj) - sm_i(ji,jj,jl1) )           ! ridge salinity = sm_i
657               sfx_bri(ji,jj) = sfx_bri(ji,jj) + sss_m(ji,jj)    * vsw(ij) * rhoic * r1_rdtice  &  ! put back sss_m into the ocean
658                  &                            - sm_i(ji,jj,jl1) * vsw(ij) * rhoic * r1_rdtice     ! and get  sm_i  from the ocean
659            ENDIF
660
661            !------------------------------------------           
662            ! 3.7 Put the snow somewhere in the ocean
663            !------------------------------------------           
664            !  Place part of the snow lost by ridging into the ocean.
665            !  Note that esrdg > 0; the ocean must cool to melt snow.
666            !  If the ocean temp = Tf already, new ice must grow.
667            !  During the next time step, thermo_rates will determine whether
668            !  the ocean cools or new ice grows.
669            zwfx_snw           = ( rhosn * vsrdg(ij) * ( 1._wp - rn_fsnowrdg )   & 
670               &                 + rhosn * vsrft(ij) * ( 1._wp - rn_fsnowrft ) ) * r1_rdtice  ! fresh water source for ocean
671
672            wfx_snw_dyn(ji,jj) = wfx_snw(ji,jj) + zwfx_snw
673            wfx_snw(ji,jj)     = wfx_snw(ji,jj) + zwfx_snw
674
675            ! SIMIP diagnostic
676            diag_dms_dyn(ji,jj) = - wfx_dyn(ji,jj)     + diag_trp_vi(ji,jj)
677            diag_dmi_dyn(ji,jj) = - wfx_snw_dyn(ji,jj) + diag_trp_vs(ji,jj)
678
679            hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ( - esrdg(ij) * ( 1._wp - rn_fsnowrdg )         & 
680               &                                - esrft(ij) * ( 1._wp - rn_fsnowrft ) ) * r1_rdtice        ! heat sink for ocean (<0, W.m-2)
681               
682            !-----------------------------------------------------------------
683            ! 3.8 Compute quantities used to apportion ice among categories
684            ! in the n2 loop below
685            !-----------------------------------------------------------------
686            dhr (ij) = 1._wp / ( hrmax(ji,jj,jl1)                    - hrmin(ji,jj,jl1)                    )
687            dhr2(ij) = 1._wp / ( hrmax(ji,jj,jl1) * hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) * hrmin(ji,jj,jl1) )
688
689
690            ! update jl1 (removing ridged/rafted area)
691            a_i  (ji,jj,  jl1) = a_i  (ji,jj,  jl1) - ardg1 (ij) - arft1 (ij)
692            v_i  (ji,jj,  jl1) = v_i  (ji,jj,  jl1) - vrdg1 (ij) - virft (ij)
693            v_s  (ji,jj,  jl1) = v_s  (ji,jj,  jl1) - vsrdg (ij) - vsrft (ij)
694            e_s  (ji,jj,1,jl1) = e_s  (ji,jj,1,jl1) - esrdg (ij) - esrft (ij)
695            smv_i(ji,jj,  jl1) = smv_i(ji,jj,  jl1) - srdg1 (ij) - smrft (ij)
696            oa_i (ji,jj,  jl1) = oa_i (ji,jj,  jl1) - oirdg1(ij) - oirft1(ij)
697
698         END DO
699
700         !--------------------------------------------------------------------
701         ! 3.9 Compute ridging ice enthalpy, remove it from ridging ice and
702         !      compute ridged ice enthalpy
703         !--------------------------------------------------------------------
704         DO jk = 1, nlay_i
705            DO ij = 1, icells
706               ji = indxi(ij) ; jj = indxj(ij)
707               ! heat content of ridged ice
708               erdg1(ij,jk) = e_i(ji,jj,jk,jl1) * afrac(ij) 
709               eirft(ij,jk) = e_i(ji,jj,jk,jl1) * afrft(ij)               
710               
711               ! enthalpy of the trapped seawater (J/m2, >0)
712               ! clem: if sst>0, then ersw <0 (is that possible?)
713               ersw(ij,jk)  = - rhoic * vsw(ij) * rcp * sst_m(ji,jj) * r1_nlay_i
714
715               ! heat flux to the ocean
716               hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ersw(ij,jk) * r1_rdtice  ! > 0 [W.m-2] ocean->ice flux
717
718               ! it is added to sea ice because the sign convention is the opposite of the sign convention for the ocean
719               erdg2(ij,jk) = erdg1(ij,jk) + ersw(ij,jk)
720
721               ! update jl1
722               e_i  (ji,jj,jk,jl1) = e_i(ji,jj,jk,jl1) - erdg1(ij,jk) - eirft(ij,jk)
723
724            END DO
725         END DO
726
727         !-------------------------------------------------------------------------------
728         ! 4) Add area, volume, and energy of new ridge to each category jl2
729         !-------------------------------------------------------------------------------
730         DO jl2  = 1, jpl 
731            ! over categories to which ridged/rafted ice is transferred
732            DO ij = 1, icells
733               ji = indxi(ij) ; jj = indxj(ij)
734
735               ! Compute the fraction of ridged ice area and volume going to thickness category jl2.
736               IF( hrmin(ji,jj,jl1) <= hi_max(jl2) .AND. hrmax(ji,jj,jl1) > hi_max(jl2-1) ) THEN
737                  hL = MAX( hrmin(ji,jj,jl1), hi_max(jl2-1) )
738                  hR = MIN( hrmax(ji,jj,jl1), hi_max(jl2)   )
739                  farea    = ( hR      - hL      ) * dhr(ij) 
740                  fvol(ij) = ( hR * hR - hL * hL ) * dhr2(ij)
741               ELSE
742                  farea    = 0._wp 
743                  fvol(ij) = 0._wp                 
744               ENDIF
745
746               ! Compute the fraction of rafted ice area and volume going to thickness category jl2
747               IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. hraft(ji,jj,jl1) >  hi_max(jl2-1) ) THEN
748                  zswitch(ij) = 1._wp
749               ELSE
750                  zswitch(ij) = 0._wp                 
751               ENDIF
752
753               a_i  (ji,jj  ,jl2) = a_i  (ji,jj  ,jl2) + ( ardg2 (ij) * farea    + arft2 (ij) * zswitch(ij) )
754               oa_i (ji,jj  ,jl2) = oa_i (ji,jj  ,jl2) + ( oirdg2(ij) * farea    + oirft2(ij) * zswitch(ij) )
755               v_i  (ji,jj  ,jl2) = v_i  (ji,jj  ,jl2) + ( vrdg2 (ij) * fvol(ij) + virft (ij) * zswitch(ij) )
756               smv_i(ji,jj  ,jl2) = smv_i(ji,jj  ,jl2) + ( srdg2 (ij) * fvol(ij) + smrft (ij) * zswitch(ij) )
757               v_s  (ji,jj  ,jl2) = v_s  (ji,jj  ,jl2) + ( vsrdg (ij) * rn_fsnowrdg * fvol(ij)  +  &
758                  &                                        vsrft (ij) * rn_fsnowrft * zswitch(ij) )
759               e_s  (ji,jj,1,jl2) = e_s  (ji,jj,1,jl2) + ( esrdg (ij) * rn_fsnowrdg * fvol(ij)  +  &
760                  &                                        esrft (ij) * rn_fsnowrft * zswitch(ij) )
761
762            END DO
763
764            ! Transfer ice energy to category jl2 by ridging
765            DO jk = 1, nlay_i
766               DO ij = 1, icells
767                  ji = indxi(ij) ; jj = indxj(ij)
768                  e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + erdg2(ij,jk) * fvol(ij) + eirft(ij,jk) * zswitch(ij)                 
769               END DO
770            END DO
771            !
772         END DO ! jl2
773         
774      END DO ! jl1 (deforming categories)
775
776      !
777      CALL wrk_dealloc( jpij,        indxi, indxj )
778      CALL wrk_dealloc( jpij,        zswitch, fvol )
779      CALL wrk_dealloc( jpij,        afrac, ardg1, ardg2, vsrdg, esrdg, dhr, dhr2 )
780      CALL wrk_dealloc( jpij,        vrdg1, vrdg2, vsw  , srdg1, srdg2, smsw, oirdg1, oirdg2 )
781      CALL wrk_dealloc( jpij,        afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 )
782      CALL wrk_dealloc( jpij,nlay_i, eirft, erdg1, erdg2, ersw )
783      !
784   END SUBROUTINE lim_itd_me_ridgeshift
785
786   SUBROUTINE lim_itd_me_icestrength( kstrngth )
787      !!----------------------------------------------------------------------
788      !!                ***  ROUTINE lim_itd_me_icestrength ***
789      !!
790      !! ** Purpose :   computes ice strength used in dynamics routines of ice thickness
791      !!
792      !! ** Method  :   Compute the strength of the ice pack, defined as the energy (J m-2)
793      !!              dissipated per unit area removed from the ice pack under compression,
794      !!              and assumed proportional to the change in potential energy caused
795      !!              by ridging. Note that only Hibler's formulation is stable and that
796      !!              ice strength has to be smoothed
797      !!
798      !! ** Inputs / Ouputs : kstrngth (what kind of ice strength we are using)
799      !!----------------------------------------------------------------------
800      INTEGER, INTENT(in) ::   kstrngth    ! = 1 for Rothrock formulation, 0 for Hibler (1979)
801      INTEGER             ::   ji,jj, jl   ! dummy loop indices
802      INTEGER             ::   ksmooth     ! smoothing the resistance to deformation
803      INTEGER             ::   numts_rm    ! number of time steps for the P smoothing
804      REAL(wp)            ::   zp, z1_3    ! local scalars
805      REAL(wp), POINTER, DIMENSION(:,:) ::   zworka   ! temporary array used here
806      !!----------------------------------------------------------------------
807
808      CALL wrk_alloc( jpi, jpj, zworka )
809
810      !------------------------------------------------------------------------------!
811      ! 1) Initialize
812      !------------------------------------------------------------------------------!
813      strength(:,:) = 0._wp
814
815      !------------------------------------------------------------------------------!
816      ! 2) Compute thickness distribution of ridging and ridged ice
817      !------------------------------------------------------------------------------!
818      CALL lim_itd_me_ridgeprep
819
820      !------------------------------------------------------------------------------!
821      ! 3) Rothrock(1975)'s method
822      !------------------------------------------------------------------------------!
823      IF( kstrngth == 1 ) THEN
824         z1_3 = 1._wp / 3._wp
825         DO jl = 1, jpl
826            DO jj= 1, jpj
827               DO ji = 1, jpi
828                  !
829                  IF( athorn(ji,jj,jl) > 0._wp ) THEN
830                     !----------------------------
831                     ! PE loss from deforming ice
832                     !----------------------------
833                     strength(ji,jj) = strength(ji,jj) - athorn(ji,jj,jl) * ht_i(ji,jj,jl) * ht_i(ji,jj,jl)
834
835                     !--------------------------
836                     ! PE gain from rafting ice
837                     !--------------------------
838                     strength(ji,jj) = strength(ji,jj) + 2._wp * araft(ji,jj,jl) * ht_i(ji,jj,jl) * ht_i(ji,jj,jl)
839
840                     !----------------------------
841                     ! PE gain from ridging ice
842                     !----------------------------
843                     strength(ji,jj) = strength(ji,jj) + aridge(ji,jj,jl) * krdg(ji,jj,jl) * z1_3 *  &
844                        &                              ( hrmax(ji,jj,jl) * hrmax(ji,jj,jl) +         &
845                        &                                hrmin(ji,jj,jl) * hrmin(ji,jj,jl) +         &
846                        &                                hrmax(ji,jj,jl) * hrmin(ji,jj,jl) ) 
847                        !!(a**3-b**3)/(a-b) = a*a+ab+b*b                     
848                  ENDIF
849                  !
850               END DO
851            END DO
852         END DO
853   
854         strength(:,:) = rn_pe_rdg * Cp * strength(:,:) / aksum(:,:)
855                         ! where Cp = (g/2)*(rhow-rhoi)*(rhoi/rhow) and rn_pe_rdg accounts for frictional dissipation
856         ksmooth = 1
857
858         !------------------------------------------------------------------------------!
859         ! 4) Hibler (1979)' method
860         !------------------------------------------------------------------------------!
861      ELSE                      ! kstrngth ne 1:  Hibler (1979) form
862         !
863         strength(:,:) = rn_pstar * vt_i(:,:) * EXP( - rn_crhg * ( 1._wp - at_i(:,:) )  )
864         !
865         ksmooth = 1
866         !
867      ENDIF                     ! kstrngth
868      !
869      !------------------------------------------------------------------------------!
870      ! 5) Impact of brine volume
871      !------------------------------------------------------------------------------!
872      ! CAN BE REMOVED
873      IF( ln_icestr_bvf ) THEN
874         DO jj = 1, jpj
875            DO ji = 1, jpi
876               strength(ji,jj) = strength(ji,jj) * exp(-5.88*SQRT(MAX(bvm_i(ji,jj),0.0)))
877            END DO
878         END DO
879      ENDIF
880      !
881      !------------------------------------------------------------------------------!
882      ! 6) Smoothing ice strength
883      !------------------------------------------------------------------------------!
884      !
885      !-------------------
886      ! Spatial smoothing
887      !-------------------
888      IF ( ksmooth == 1 ) THEN
889
890         CALL lbc_lnk( strength, 'T', 1. )
891
892         DO jj = 2, jpjm1
893            DO ji = 2, jpim1
894               IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp) THEN
895                  zworka(ji,jj) = ( 4.0 * strength(ji,jj)              &
896                     &                  + strength(ji-1,jj) * tmask(ji-1,jj,1) + strength(ji+1,jj) * tmask(ji+1,jj,1) & 
897                     &                  + strength(ji,jj-1) * tmask(ji,jj-1,1) + strength(ji,jj+1) * tmask(ji,jj+1,1) &
898                     &            ) / ( 4.0 + tmask(ji-1,jj,1) + tmask(ji+1,jj,1) + tmask(ji,jj-1,1) + tmask(ji,jj+1,1) )
899               ELSE
900                  zworka(ji,jj) = 0._wp
901               ENDIF
902            END DO
903         END DO
904
905         DO jj = 2, jpjm1
906            DO ji = 2, jpim1
907               strength(ji,jj) = zworka(ji,jj)
908            END DO
909         END DO
910         CALL lbc_lnk( strength, 'T', 1. )
911
912      ENDIF
913
914      !--------------------
915      ! Temporal smoothing
916      !--------------------
917      IF ( numit == nit000 + nn_fsbc - 1 ) THEN
918         strp1(:,:) = 0.0           
919         strp2(:,:) = 0.0           
920      ENDIF
921
922      IF ( ksmooth == 2 ) THEN
923
924         CALL lbc_lnk( strength, 'T', 1. )
925
926         DO jj = 1, jpj - 1
927            DO ji = 1, jpi - 1
928               IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp) THEN
929                  numts_rm = 1 ! number of time steps for the running mean
930                  IF ( strp1(ji,jj) > 0._wp ) numts_rm = numts_rm + 1
931                  IF ( strp2(ji,jj) > 0._wp ) numts_rm = numts_rm + 1
932                  zp = ( strength(ji,jj) + strp1(ji,jj) + strp2(ji,jj) ) / numts_rm
933                  strp2(ji,jj) = strp1(ji,jj)
934                  strp1(ji,jj) = strength(ji,jj)
935                  strength(ji,jj) = zp
936
937               ENDIF
938            END DO
939         END DO
940
941      ENDIF ! ksmooth
942
943      CALL lbc_lnk( strength, 'T', 1. )      ! Boundary conditions
944
945      CALL wrk_dealloc( jpi, jpj, zworka )
946      !
947   END SUBROUTINE lim_itd_me_icestrength
948
949   SUBROUTINE lim_itd_me_init
950      !!-------------------------------------------------------------------
951      !!                   ***  ROUTINE lim_itd_me_init ***
952      !!
953      !! ** Purpose :   Physical constants and parameters linked
954      !!                to the mechanical ice redistribution
955      !!
956      !! ** Method  :   Read the namiceitdme namelist
957      !!                and check the parameters values
958      !!                called at the first timestep (nit000)
959      !!
960      !! ** input   :   Namelist namiceitdme
961      !!-------------------------------------------------------------------
962      INTEGER :: ios                 ! Local integer output status for namelist read
963      NAMELIST/namiceitdme/ rn_cs, rn_fsnowrdg, rn_fsnowrft,              & 
964        &                   rn_gstar, rn_astar, rn_hstar, ln_rafting, rn_hraft, rn_craft, rn_por_rdg, &
965        &                   nn_partfun
966      !!-------------------------------------------------------------------
967      !
968      REWIND( numnam_ice_ref )              ! Namelist namicetdme in reference namelist : Ice mechanical ice redistribution
969      READ  ( numnam_ice_ref, namiceitdme, IOSTAT = ios, ERR = 901)
970901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitdme in reference namelist', lwp )
971
972      REWIND( numnam_ice_cfg )              ! Namelist namiceitdme in configuration namelist : Ice mechanical ice redistribution
973      READ  ( numnam_ice_cfg, namiceitdme, IOSTAT = ios, ERR = 902 )
974902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitdme in configuration namelist', lwp )
975      IF(lwm) WRITE ( numoni, namiceitdme )
976      !
977      IF (lwp) THEN                          ! control print
978         WRITE(numout,*)
979         WRITE(numout,*)' lim_itd_me_init : ice parameters for mechanical ice redistribution '
980         WRITE(numout,*)' ~~~~~~~~~~~~~~~'
981         WRITE(numout,*)'   Fraction of shear energy contributing to ridging        rn_cs       = ', rn_cs 
982         WRITE(numout,*)'   Fraction of snow volume conserved during ridging        rn_fsnowrdg = ', rn_fsnowrdg 
983         WRITE(numout,*)'   Fraction of snow volume conserved during ridging        rn_fsnowrft = ', rn_fsnowrft 
984         WRITE(numout,*)'   Fraction of total ice coverage contributing to ridging  rn_gstar    = ', rn_gstar
985         WRITE(numout,*)'   Equivalent to G* for an exponential part function       rn_astar    = ', rn_astar
986         WRITE(numout,*)'   Quantity playing a role in max ridged ice thickness     rn_hstar    = ', rn_hstar
987         WRITE(numout,*)'   Rafting of ice sheets or not                            ln_rafting  = ', ln_rafting
988         WRITE(numout,*)'   Parmeter thickness (threshold between ridge-raft)       rn_hraft    = ', rn_hraft
989         WRITE(numout,*)'   Rafting hyperbolic tangent coefficient                  rn_craft    = ', rn_craft 
990         WRITE(numout,*)'   Initial porosity of ridges                              rn_por_rdg  = ', rn_por_rdg
991         WRITE(numout,*)'   Switch for part. function (0) linear (1) exponential    nn_partfun  = ', nn_partfun
992      ENDIF
993      !
994   END SUBROUTINE lim_itd_me_init
995
996#else
997   !!----------------------------------------------------------------------
998   !!   Default option         Empty module          NO LIM-3 sea-ice model
999   !!----------------------------------------------------------------------
1000CONTAINS
1001   SUBROUTINE lim_itd_me           ! Empty routines
1002   END SUBROUTINE lim_itd_me
1003   SUBROUTINE lim_itd_me_icestrength
1004   END SUBROUTINE lim_itd_me_icestrength
1005   SUBROUTINE lim_itd_me_init
1006   END SUBROUTINE lim_itd_me_init
1007#endif
1008   !!======================================================================
1009END MODULE limitd_me
Note: See TracBrowser for help on using the repository browser.