source: branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 @ 8142

Last change on this file since 8142 was 8142, checked in by vancop, 3 years ago

Melt pond interfaces practically operational

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