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

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

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