New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
limitd_me.F90 in branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

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

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

a few more melt pond corrections

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