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

source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 @ 8319

Last change on this file since 8319 was 8319, checked in by clem, 7 years ago

STEP2 (4): remove obsolete features (change time step for the ice => remove numit, nstart, nlast, nitrun and create kt_ice which follows kt)

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