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.
icerdgrft.F90 in branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerdgrft.F90 @ 8426

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

last routine names to be changed

File size: 48.9 KB
Line 
1MODULE icerdgrft
2   !!======================================================================
3   !!                       ***  MODULE icerdgrft ***
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 ice1D            ! LIM thermodynamics
19   USE ice              ! LIM variables
20   USE icevar           ! LIM
21   USE icectl           ! control prints
22   !
23   USE lbclnk           ! lateral boundary condition - MPP exchanges
24   USE lib_mpp          ! MPP library
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
30   IMPLICIT NONE
31   PRIVATE
32
33   PUBLIC   ice_rdgrft               ! called by ice_stp
34   PUBLIC   ice_rdgrft_icestrength
35   PUBLIC   ice_rdgrft_init
36   PUBLIC   ice_rdgrft_alloc        ! called by ice_init
37
38   !-----------------------------------------------------------------------
39   ! Variables shared among ridging subroutines
40   !-----------------------------------------------------------------------
41   REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   asum     ! sum of total ice and open water area
42   REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   aksum    ! ratio of area removed to area ridged
43   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   athorn   ! participation function; fraction of ridging/closing associated w/ category n
44   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   hrmin    ! minimum ridge thickness
45   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   hrmax    ! maximum ridge thickness
46   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   hraft    ! thickness of rafted ice
47   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   krdg     ! thickness of ridging ice / mean ridge thickness
48   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   aridge   ! participating ice ridging
49   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   araft    ! participating ice rafting
50
51   REAL(wp), PARAMETER ::   krdgmin = 1.1_wp    ! min ridge thickness multiplier
52   REAL(wp), PARAMETER ::   kraft   = 0.5_wp    ! rafting multipliyer
53
54   REAL(wp) ::   Cp                             !
55   !
56   !
57   !!----------------------------------------------------------------------
58   !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010)
59   !! $Id: icerdgrft.F90 8378 2017-07-26 13:55:59Z clem $
60   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
61   !!----------------------------------------------------------------------
62CONTAINS
63
64   INTEGER FUNCTION ice_rdgrft_alloc()
65      !!---------------------------------------------------------------------!
66      !!                ***  ROUTINE ice_rdgrft_alloc ***
67      !!---------------------------------------------------------------------!
68      ALLOCATE(                                                                      &
69         !* Variables shared among ridging subroutines
70         &      asum (jpi,jpj)     , athorn(jpi,jpj,0:jpl) , aksum (jpi,jpj)     ,   &
71         &      hrmin(jpi,jpj,jpl) , hraft(jpi,jpj,jpl)    , aridge(jpi,jpj,jpl) ,   &
72         &      hrmax(jpi,jpj,jpl) , krdg (jpi,jpj,jpl)    , araft (jpi,jpj,jpl) , STAT=ice_rdgrft_alloc )
73         !
74      IF( ice_rdgrft_alloc /= 0 )   CALL ctl_warn( 'ice_rdgrft_alloc: failed to allocate arrays' )
75      !
76   END FUNCTION ice_rdgrft_alloc
77
78
79   SUBROUTINE ice_rdgrft( kt )
80      !!---------------------------------------------------------------------!
81      !!                ***  ROUTINE ice_rdgrft ***
82      !!
83      !! ** Purpose :   computes the mechanical redistribution of ice thickness
84      !!
85      !! ** Method  :   Steps :
86      !!       1) Thickness categories boundaries, ice / o.w. concentrations
87      !!          Ridge preparation
88      !!       2) Dynamical inputs (closing rate, divu_adv, opning)
89      !!       3) Ridging iteration
90      !!       4) Ridging diagnostics
91      !!       5) Heat, salt and freshwater fluxes
92      !!       6) Compute increments of tate variables and come back to old values
93      !!
94      !! References :   Flato, G. M., and W. D. Hibler III, 1995, JGR, 100, 18,611-18,626.
95      !!                Hibler, W. D. III, 1980, MWR, 108, 1943-1973, 1980.
96      !!                Rothrock, D. A., 1975: JGR, 80, 4514-4519.
97      !!                Thorndike et al., 1975, JGR, 80, 4501-4513.
98      !!                Bitz et al., JGR, 2001
99      !!                Amundrud and Melling, JGR 2005
100      !!                Babko et al., JGR 2002
101      !!
102      !!     This routine is based on CICE code and authors William H. Lipscomb,
103      !!  and Elizabeth C. Hunke, LANL are gratefully acknowledged
104      !!--------------------------------------------------------------------!
105      INTEGER, INTENT(in) ::   kt     ! number of iteration
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), DIMENSION(jpi,jpj)   ::   closing_net     ! net rate at which area is removed    (1/s)
113                                                               ! (ridging ice area - area of new ridges) / dt
114      REAL(wp), DIMENSION(jpi,jpj)   ::   divu_adv        ! divu as implied by transport scheme  (1/s)
115      REAL(wp), DIMENSION(jpi,jpj)   ::   opning          ! rate of opening due to divergence/shear
116      REAL(wp), DIMENSION(jpi,jpj)   ::   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('icerdgrft')
123
124      IF( kt == nit000 .AND. lwp ) THEN
125         WRITE(numout,*)
126         WRITE(numout,*)'icerdgrft'
127         WRITE(numout,*)'~~~~~~~~~'
128      ENDIF
129
130      ! conservation test
131      IF( ln_limdiachk ) CALL ice_cons_hsm(0, 'icerdgrft', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
132
133      !-----------------------------------------------------------------------------!
134      ! 1) Thickness categories boundaries, ice / o.w. concentrations, init_ons
135      !-----------------------------------------------------------------------------!
136      Cp = 0.5 * grav * (rau0-rhoic) * rhoic * r1_rau0             ! proport const for PE
137      !
138      CALL ice_rdgrft_ridgeprep                                    ! prepare ridging
139      !
140
141      DO jj = 1, jpj                                     ! Initialize arrays.
142         DO ji = 1, jpi
143
144            !-----------------------------------------------------------------------------!
145            ! 2) Dynamical inputs (closing rate, divu_adv, opning)
146            !-----------------------------------------------------------------------------!
147            !
148            ! 2.1 closing_net
149            !-----------------
150            ! Compute the net rate of closing due to convergence
151            ! and shear, based on Flato and Hibler (1995).
152            !
153            ! The energy dissipation rate is equal to the net closing rate
154            ! times the ice strength.
155            !
156            ! NOTE: The NET closing rate is equal to the rate that open water
157            !  area is removed, plus the rate at which ice area is removed by
158            !  ridging, minus the rate at which area is added in new ridges.
159            !  The GROSS closing rate is equal to the first two terms (open
160            !  water closing and thin ice ridging) without the third term
161            !  (thick, newly ridged ice).
162
163            closing_net(ji,jj) = rn_cs * 0.5 * ( delta_i(ji,jj) - ABS( divu_i(ji,jj) ) ) - MIN( divu_i(ji,jj), 0._wp )
164
165            ! 2.2 divu_adv
166            !--------------
167            ! Compute divu_adv, the divergence rate given by the transport/
168            ! advection scheme, which may not be equal to divu as computed
169            ! from the velocity field.
170            !
171            ! If divu_adv < 0, make sure the closing rate is large enough
172            ! to give asum = 1.0 after ridging.
173           
174            divu_adv(ji,jj) = ( 1._wp - asum(ji,jj) ) * r1_rdtice  ! asum found in ridgeprep
175
176            IF( divu_adv(ji,jj) < 0._wp )   closing_net(ji,jj) = MAX( closing_net(ji,jj), -divu_adv(ji,jj) )
177
178            ! 2.3 opning
179            !------------
180            ! Compute the (non-negative) opening rate that will give asum = 1.0 after ridging.
181            opning(ji,jj) = closing_net(ji,jj) + divu_adv(ji,jj)
182         END DO
183      END DO
184
185      !-----------------------------------------------------------------------------!
186      ! 3) Ridging iteration
187      !-----------------------------------------------------------------------------!
188      niter           = 1                 ! iteration counter
189      iterate_ridging = 1
190
191      DO WHILE ( iterate_ridging > 0 .AND. niter < nitermax )
192
193         ! 3.2 closing_gross
194         !-----------------------------------------------------------------------------!
195         ! Based on the ITD of ridging and ridged ice, convert the net
196         !  closing rate to a gross closing rate. 
197         ! NOTE: 0 < aksum <= 1
198         closing_gross(:,:) = closing_net(:,:) / aksum(:,:)
199
200         ! correction to closing rate and opening if closing rate is excessive
201         !---------------------------------------------------------------------
202         ! Reduce the closing rate if more than 100% of the open water
203         ! would be removed.  Reduce the opening rate proportionately.
204         DO jj = 1, jpj
205            DO ji = 1, jpi
206               za   = ( opning(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) ) * rdt_ice
207               IF    ( za < 0._wp .AND. za > - ato_i(ji,jj) ) THEN                  ! would lead to negative ato_i
208                  zfac          = - ato_i(ji,jj) / za
209                  opning(ji,jj) = athorn(ji,jj,0) * closing_gross(ji,jj) - ato_i(ji,jj) * r1_rdtice 
210               ELSEIF( za > 0._wp .AND. za > ( asum(ji,jj) - ato_i(ji,jj) ) ) THEN  ! would lead to ato_i > asum
211                  zfac          = ( asum(ji,jj) - ato_i(ji,jj) ) / za
212                  opning(ji,jj) = athorn(ji,jj,0) * closing_gross(ji,jj) + ( asum(ji,jj) - ato_i(ji,jj) ) * r1_rdtice 
213               ENDIF
214            END DO
215         END DO
216
217         ! correction to closing rate / opening if excessive ice removal
218         !---------------------------------------------------------------
219         ! Reduce the closing rate if more than 100% of any ice category
220         ! would be removed.  Reduce the opening rate proportionately.
221         DO jl = 1, jpl
222            DO jj = 1, jpj
223               DO ji = 1, jpi
224                  za = athorn(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice
225                  IF( za  >  a_i(ji,jj,jl) ) THEN
226                     zfac = a_i(ji,jj,jl) / za
227                     closing_gross(ji,jj) = closing_gross(ji,jj) * zfac
228                  ENDIF
229               END DO
230            END DO
231         END DO
232
233         ! 3.3 Redistribute area, volume, and energy.
234         !-----------------------------------------------------------------------------!
235
236         CALL ice_rdgrft_ridgeshift( opning, closing_gross )
237
238         
239         ! 3.4 Compute total area of ice plus open water after ridging.
240         !-----------------------------------------------------------------------------!
241         ! This is in general not equal to one because of divergence during transport
242         asum(:,:) = ato_i(:,:) + SUM( a_i, dim=3 )
243
244         ! 3.5 Do we keep on iterating ???
245         !-----------------------------------------------------------------------------!
246         ! Check whether asum = 1.  If not (because the closing and opening
247         ! rates were reduced above), ridge again with new rates.
248
249         iterate_ridging = 0
250         DO jj = 1, jpj
251            DO ji = 1, jpi
252               IF( ABS( asum(ji,jj) - 1._wp ) < epsi10 ) THEN
253                  closing_net(ji,jj) = 0._wp
254                  opning     (ji,jj) = 0._wp
255                  ato_i      (ji,jj) = MAX( 0._wp, 1._wp - SUM( a_i(ji,jj,:) ) )
256               ELSE
257                  iterate_ridging    = 1
258                  divu_adv   (ji,jj) = ( 1._wp - asum(ji,jj) ) * r1_rdtice
259                  closing_net(ji,jj) = MAX( 0._wp, -divu_adv(ji,jj) )
260                  opning     (ji,jj) = MAX( 0._wp,  divu_adv(ji,jj) )
261               ENDIF
262            END DO
263         END DO
264
265         IF( lk_mpp )   CALL mpp_max( iterate_ridging )
266
267         ! Repeat if necessary.
268         ! NOTE: If strength smoothing is turned on, the ridging must be
269         ! iterated globally because of the boundary update in the
270         ! smoothing.
271
272         niter = niter + 1
273
274         IF( iterate_ridging == 1 ) THEN
275            CALL ice_rdgrft_ridgeprep
276            IF( niter  >  nitermax ) THEN
277               WRITE(numout,*) ' ALERTE : non-converging ridging scheme '
278               WRITE(numout,*) ' niter, iterate_ridging ', niter, iterate_ridging
279            ENDIF
280         ENDIF
281
282      END DO !! on the do while over iter
283
284      CALL ice_var_agg( 1 ) 
285
286      !-----------------------------------------------------------------------------!
287      ! control prints
288      !-----------------------------------------------------------------------------!
289      ! conservation test
290      IF( ln_limdiachk ) CALL ice_cons_hsm(1, 'icerdgrft', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
291
292      ! control prints
293      IF( ln_ctl )       CALL ice_prt3D( 'icerdgrft' )
294      !
295      IF( nn_timing == 1 )  CALL timing_stop('icerdgrft')
296   END SUBROUTINE ice_rdgrft
297
298   SUBROUTINE ice_rdgrft_ridgeprep
299      !!---------------------------------------------------------------------!
300      !!                ***  ROUTINE ice_rdgrft_ridgeprep ***
301      !!
302      !! ** Purpose :   preparation for ridging and strength calculations
303      !!
304      !! ** Method  :   Compute the thickness distribution of the ice and open water
305      !!              participating in ridging and of the resulting ridges.
306      !!---------------------------------------------------------------------!
307      INTEGER  ::   ji,jj, jl    ! dummy loop indices
308      REAL(wp) ::   Gstari, astari, hrmean, zdummy   ! local scalar
309      REAL(wp), DIMENSION(jpi,jpj,-1:jpl) ::   Gsum      ! Gsum(n) = sum of areas in categories 0 to n
310      !------------------------------------------------------------------------------!
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 ice_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      !
480   END SUBROUTINE ice_rdgrft_ridgeprep
481
482
483   SUBROUTINE ice_rdgrft_ridgeshift( opning, closing_gross )
484      !!----------------------------------------------------------------------
485      !!                ***  ROUTINE ice_rdgrft_icestrength ***
486      !!
487      !! ** Purpose :   shift ridging ice among thickness categories of ice thickness
488      !!
489      !! ** Method  :   Remove area, volume, and energy from each ridging category
490      !!              and add to thicker ice categories.
491      !!----------------------------------------------------------------------
492      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   opning         ! rate of opening due to divergence/shear
493      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   closing_gross  ! rate at which area removed, excluding area of new ridges
494      !
495      CHARACTER (len=80) ::   fieldid   ! field identifier
496      !
497      INTEGER ::   ji, jj, jl, jl1, jl2, jk   ! dummy loop indices
498      INTEGER ::   ij                ! horizontal index, combines i and j loops
499      INTEGER ::   icells            ! number of cells with a_i > puny
500      REAL(wp) ::   hL, hR, farea    ! left and right limits of integration
501      REAL(wp) ::   zwfx_snw         ! snow mass flux increment
502
503      INTEGER , DIMENSION(jpij) ::   indxi, indxj   ! compressed indices
504      REAL(wp), DIMENSION(jpij) ::   zswitch, fvol   ! new ridge volume going to n2
505
506      REAL(wp), DIMENSION(jpij) ::   afrac            ! fraction of category area ridged
507      REAL(wp), DIMENSION(jpij) ::   ardg1 , ardg2    ! area of ice ridged & new ridges
508      REAL(wp), DIMENSION(jpij) ::   vsrdg , esrdg    ! snow volume & energy of ridging ice
509      ! MV MP 2016
510      REAL(wp), DIMENSION(jpij) ::   vprdg            ! pond volume of ridging ice
511      REAL(wp), DIMENSION(jpij) ::   aprdg1           ! pond area of ridging ice
512      REAL(wp), DIMENSION(jpij) ::   aprdg2           ! pond area of ridging ice
513      ! END MV MP 2016
514      REAL(wp), DIMENSION(jpij) ::   dhr   , dhr2     ! hrmax - hrmin  &  hrmax^2 - hrmin^2
515
516      REAL(wp), DIMENSION(jpij) ::   vrdg1   ! volume of ice ridged
517      REAL(wp), DIMENSION(jpij) ::   vrdg2   ! volume of new ridges
518      REAL(wp), DIMENSION(jpij) ::   vsw     ! volume of seawater trapped into ridges
519      REAL(wp), DIMENSION(jpij) ::   srdg1   ! sal*volume of ice ridged
520      REAL(wp), DIMENSION(jpij) ::   srdg2   ! sal*volume of new ridges
521      REAL(wp), DIMENSION(jpij) ::   smsw    ! sal*volume of water trapped into ridges
522      REAL(wp), DIMENSION(jpij) ::   oirdg1, oirdg2   ! ice age of ice ridged
523
524      REAL(wp), DIMENSION(jpij) ::   afrft            ! fraction of category area rafted
525      REAL(wp), DIMENSION(jpij) ::   arft1 , arft2    ! area of ice rafted and new rafted zone
526      REAL(wp), DIMENSION(jpij) ::   virft , vsrft    ! ice & snow volume of rafting ice
527      ! MV MP 2016
528      REAL(wp), DIMENSION(jpij) ::   vprft            ! pond volume of rafting ice
529      REAL(wp), DIMENSION(jpij) ::   aprft1           ! pond area of rafted ice
530      REAL(wp), DIMENSION(jpij) ::   aprft2           ! pond area of new rafted ice
531      ! END MV MP 2016
532      REAL(wp), DIMENSION(jpij) ::   esrft , smrft    ! snow energy & salinity of rafting ice
533      REAL(wp), DIMENSION(jpij) ::   oirft1, oirft2   ! ice age of ice rafted
534
535      REAL(wp), DIMENSION(jpij,nlay_i) ::   eirft      ! ice energy of rafting ice
536      REAL(wp), DIMENSION(jpij,nlay_i) ::   erdg1      ! enth*volume of ice ridged
537      REAL(wp), DIMENSION(jpij,nlay_i) ::   erdg2      ! enth*volume of new ridges
538      REAL(wp), DIMENSION(jpij,nlay_i) ::   ersw       ! enth of water trapped into ridges
539      !!----------------------------------------------------------------------
540
541      !-------------------------------------------------------------------------------
542      ! 1) Compute change in open water area due to closing and opening.
543      !-------------------------------------------------------------------------------
544      DO jj = 1, jpj
545         DO ji = 1, jpi
546            ato_i(ji,jj) = MAX( 0._wp, ato_i(ji,jj) +  &
547               &                     ( opning(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) ) * rdt_ice )
548         END DO
549      END DO
550
551      !-----------------------------------------------------------------
552      ! 3) Pump everything from ice which is being ridged / rafted
553      !-----------------------------------------------------------------
554      ! Compute the area, volume, and energy of ice ridging in each
555      ! category, along with the area of the resulting ridge.
556
557      DO jl1 = 1, jpl !jl1 describes the ridging category
558
559         !------------------------------------------------
560         ! 3.1) Identify grid cells with nonzero ridging
561         !------------------------------------------------
562         icells = 0
563         DO jj = 1, jpj
564            DO ji = 1, jpi
565               IF( athorn(ji,jj,jl1) > 0._wp .AND. closing_gross(ji,jj) > 0._wp ) THEN
566                  icells = icells + 1
567                  indxi(icells) = ji
568                  indxj(icells) = jj
569               ENDIF
570            END DO
571         END DO
572
573         DO ij = 1, icells
574            ji = indxi(ij) ; jj = indxj(ij)
575
576            !--------------------------------------------------------------------
577            ! 3.2) Compute area of ridging ice (ardg1) and of new ridge (ardg2)
578            !--------------------------------------------------------------------
579            ardg1(ij) = aridge(ji,jj,jl1) * closing_gross(ji,jj) * rdt_ice
580            arft1(ij) = araft (ji,jj,jl1) * closing_gross(ji,jj) * rdt_ice
581
582            !---------------------------------------------------------------
583            ! 3.3) Compute ridging /rafting fractions, make sure afrac <=1
584            !---------------------------------------------------------------
585            afrac(ij) = ardg1(ij) / a_i(ji,jj,jl1) !ridging
586            afrft(ij) = arft1(ij) / a_i(ji,jj,jl1) !rafting
587            ardg2(ij) = ardg1(ij) * krdg(ji,jj,jl1)
588            arft2(ij) = arft1(ij) * kraft
589
590            !--------------------------------------------------------------------------
591            ! 3.4) Substract area, volume, and energy from ridging
592            !     / rafting category n1.
593            !--------------------------------------------------------------------------
594            vrdg1(ij) = v_i(ji,jj,jl1) * afrac(ij)
595            vrdg2(ij) = vrdg1(ij) * ( 1. + rn_por_rdg )
596            vsw  (ij) = vrdg1(ij) * rn_por_rdg
597
598            vsrdg (ij) = v_s  (ji,jj,  jl1) * afrac(ij)
599            esrdg (ij) = e_s  (ji,jj,1,jl1) * afrac(ij)
600            !MV MP 2016
601            IF ( nn_pnd_scheme > 0 ) THEN
602               aprdg1(ij) = a_ip(ji,jj, jl1) * afrac(ij)
603               aprdg2(ij) = a_ip(ji,jj, jl1) * afrac(ij) * krdg(ji,jj,jl1)
604               vprdg(ij)  = v_ip(ji,jj, jl1) * afrac(ij)
605            ENDIF
606            ! END MV MP 2016
607            srdg1 (ij) = smv_i(ji,jj,  jl1) * afrac(ij)
608            oirdg1(ij) = oa_i (ji,jj,  jl1) * afrac(ij)
609            oirdg2(ij) = oa_i (ji,jj,  jl1) * afrac(ij) * krdg(ji,jj,jl1) 
610
611            ! rafting volumes, heat contents ...
612            virft (ij) = v_i  (ji,jj,  jl1) * afrft(ij)
613            vsrft (ij) = v_s  (ji,jj,  jl1) * afrft(ij)
614            !MV MP 2016
615            IF ( nn_pnd_scheme > 0 ) THEN
616               aprft1(ij) = a_ip (ji,jj,  jl1) * afrft(ij)
617               aprft2(ij) = a_ip (ji,jj,  jl1) * afrft(ij) * kraft
618               vprft(ij)  = v_ip(ji,jj,jl1)    * afrft(ij)
619            ENDIF
620            ! END MV MP 2016
621            srdg1 (ij) = smv_i(ji,jj,  jl1) * afrac(ij)
622            esrft (ij) = e_s  (ji,jj,1,jl1) * afrft(ij)
623            smrft (ij) = smv_i(ji,jj,  jl1) * afrft(ij) 
624            oirft1(ij) = oa_i (ji,jj,  jl1) * afrft(ij) 
625            oirft2(ij) = oa_i (ji,jj,  jl1) * afrft(ij) * kraft 
626
627            !-----------------------------------------------------------------
628            ! 3.5) Compute properties of new ridges
629            !-----------------------------------------------------------------
630            smsw(ij)  = vsw(ij) * sss_m(ji,jj)                   ! salt content of seawater frozen in voids
631            srdg2(ij) = srdg1(ij) + smsw(ij)                     ! salt content of new ridge
632           
633            sfx_dyn(ji,jj) = sfx_dyn(ji,jj) - smsw(ij) * rhoic * r1_rdtice
634            wfx_dyn(ji,jj) = wfx_dyn(ji,jj) - vsw (ij) * rhoic * r1_rdtice   ! increase in ice volume due to seawater frozen in voids
635
636            ! virtual salt flux to keep salinity constant
637            IF( nn_icesal == 1 .OR. nn_icesal == 3 )  THEN
638               srdg2(ij)      = srdg2(ij) - vsw(ij) * ( sss_m(ji,jj) - sm_i(ji,jj,jl1) )           ! ridge salinity = sm_i
639               sfx_bri(ji,jj) = sfx_bri(ji,jj) + sss_m(ji,jj)    * vsw(ij) * rhoic * r1_rdtice  &  ! put back sss_m into the ocean
640                  &                            - sm_i(ji,jj,jl1) * vsw(ij) * rhoic * r1_rdtice     ! and get  sm_i  from the ocean
641            ENDIF
642
643            !------------------------------------------           
644            ! 3.7 Put the snow somewhere in the ocean
645            !------------------------------------------           
646            !  Place part of the snow lost by ridging into the ocean.
647            !  Note that esrdg > 0; the ocean must cool to melt snow.
648            !  If the ocean temp = Tf already, new ice must grow.
649            !  During the next time step, thermo_rates will determine whether
650            !  the ocean cools or new ice grows.
651            zwfx_snw           = ( rhosn * vsrdg(ij) * ( 1._wp - rn_fsnowrdg )   & 
652               &                 + rhosn * vsrft(ij) * ( 1._wp - rn_fsnowrft ) ) * r1_rdtice  ! fresh water source for ocean
653
654            wfx_snw_dyn(ji,jj)  =   wfx_snw_dyn(ji,jj) + zwfx_snw
655
656            hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ( - esrdg(ij) * ( 1._wp - rn_fsnowrdg )         & 
657               &                                - esrft(ij) * ( 1._wp - rn_fsnowrft ) ) * r1_rdtice        ! heat sink for ocean (<0, W.m-2)
658
659            ! MV MP 2016
660            !------------------------------------------           
661            ! 3.X Put the melt pond water in the ocean
662            !------------------------------------------           
663            !  Place part of the melt pond volume into the ocean.
664            IF ( ( nn_pnd_scheme > 0 ) .AND. ln_pnd_fw ) THEN
665               wfx_pnd(ji,jj) = wfx_pnd(ji,jj) + ( rhofw * vprdg(ij) * ( 1._wp - rn_fpondrdg )   & 
666               &                                 + rhofw * vprft(ij) * ( 1._wp - rn_fpondrft ) ) * r1_rdtice  ! fresh water source for ocean
667            ENDIF
668            ! END MV MP 2016
669
670            !-----------------------------------------------------------------
671            ! 3.8 Compute quantities used to apportion ice among categories
672            ! in the n2 loop below
673            !-----------------------------------------------------------------
674            dhr (ij) = 1._wp / ( hrmax(ji,jj,jl1)                    - hrmin(ji,jj,jl1)                    )
675            dhr2(ij) = 1._wp / ( hrmax(ji,jj,jl1) * hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) * hrmin(ji,jj,jl1) )
676
677
678            ! update jl1 (removing ridged/rafted area)
679            a_i  (ji,jj,  jl1) = a_i  (ji,jj,  jl1) - ardg1 (ij) - arft1 (ij)
680            v_i  (ji,jj,  jl1) = v_i  (ji,jj,  jl1) - vrdg1 (ij) - virft (ij)
681            v_s  (ji,jj,  jl1) = v_s  (ji,jj,  jl1) - vsrdg (ij) - vsrft (ij)
682            e_s  (ji,jj,1,jl1) = e_s  (ji,jj,1,jl1) - esrdg (ij) - esrft (ij)
683            smv_i(ji,jj,  jl1) = smv_i(ji,jj,  jl1) - srdg1 (ij) - smrft (ij)
684            oa_i (ji,jj,  jl1) = oa_i (ji,jj,  jl1) - oirdg1(ij) - oirft1(ij)
685
686            ! MV MP 2016
687            IF ( nn_pnd_scheme > 0 ) THEN
688               v_ip (ji,jj,jl1) = v_ip (ji,jj,jl1) - vprdg (ij) - vprft (ij)
689               a_ip (ji,jj,jl1) = a_ip (ji,jj,jl1) - aprdg1(ij) - aprft1(ij)
690            ENDIF
691            ! END MV MP 2016
692
693         END DO
694
695         !--------------------------------------------------------------------
696         ! 3.9 Compute ridging ice enthalpy, remove it from ridging ice and
697         !      compute ridged ice enthalpy
698         !--------------------------------------------------------------------
699         DO jk = 1, nlay_i
700            DO ij = 1, icells
701               ji = indxi(ij) ; jj = indxj(ij)
702               ! heat content of ridged ice
703               erdg1(ij,jk) = e_i(ji,jj,jk,jl1) * afrac(ij) 
704               eirft(ij,jk) = e_i(ji,jj,jk,jl1) * afrft(ij)               
705               
706               ! enthalpy of the trapped seawater (J/m2, >0)
707               ! clem: if sst>0, then ersw <0 (is that possible?)
708               ersw(ij,jk)  = - rhoic * vsw(ij) * rcp * sst_m(ji,jj) * r1_nlay_i
709
710               ! heat flux to the ocean
711               hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ersw(ij,jk) * r1_rdtice  ! > 0 [W.m-2] ocean->ice flux
712
713               ! it is added to sea ice because the sign convention is the opposite of the sign convention for the ocean
714               erdg2(ij,jk) = erdg1(ij,jk) + ersw(ij,jk)
715
716               ! update jl1
717               e_i  (ji,jj,jk,jl1) = e_i(ji,jj,jk,jl1) - erdg1(ij,jk) - eirft(ij,jk)
718
719            END DO
720         END DO
721
722         !-------------------------------------------------------------------------------
723         ! 4) Add area, volume, and energy of new ridge to each category jl2
724         !-------------------------------------------------------------------------------
725         DO jl2  = 1, jpl 
726            ! over categories to which ridged/rafted ice is transferred
727            DO ij = 1, icells
728               ji = indxi(ij) ; jj = indxj(ij)
729
730               ! Compute the fraction of ridged ice area and volume going to thickness category jl2.
731               IF( hrmin(ji,jj,jl1) <= hi_max(jl2) .AND. hrmax(ji,jj,jl1) > hi_max(jl2-1) ) THEN
732                  hL = MAX( hrmin(ji,jj,jl1), hi_max(jl2-1) )
733                  hR = MIN( hrmax(ji,jj,jl1), hi_max(jl2)   )
734                  farea    = ( hR      - hL      ) * dhr(ij) 
735                  fvol(ij) = ( hR * hR - hL * hL ) * dhr2(ij)
736               ELSE
737                  farea    = 0._wp 
738                  fvol(ij) = 0._wp                 
739               ENDIF
740
741               ! Compute the fraction of rafted ice area and volume going to thickness category jl2
742               IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. hraft(ji,jj,jl1) >  hi_max(jl2-1) ) THEN
743                  zswitch(ij) = 1._wp
744               ELSE
745                  zswitch(ij) = 0._wp                 
746               ENDIF
747
748               a_i  (ji,jj  ,jl2) = a_i  (ji,jj  ,jl2) + ( ardg2 (ij) * farea    + arft2 (ij) * zswitch(ij) )
749               oa_i (ji,jj  ,jl2) = oa_i (ji,jj  ,jl2) + ( oirdg2(ij) * farea    + oirft2(ij) * zswitch(ij) )
750               v_i  (ji,jj  ,jl2) = v_i  (ji,jj  ,jl2) + ( vrdg2 (ij) * fvol(ij) + virft (ij) * zswitch(ij) )
751               smv_i(ji,jj  ,jl2) = smv_i(ji,jj  ,jl2) + ( srdg2 (ij) * fvol(ij) + smrft (ij) * zswitch(ij) )
752               v_s  (ji,jj  ,jl2) = v_s  (ji,jj  ,jl2) + ( vsrdg (ij) * rn_fsnowrdg * fvol(ij)  +  &
753                  &                                        vsrft (ij) * rn_fsnowrft * zswitch(ij) )
754               e_s  (ji,jj,1,jl2) = e_s  (ji,jj,1,jl2) + ( esrdg (ij) * rn_fsnowrdg * fvol(ij)  +  &
755                  &                                        esrft (ij) * rn_fsnowrft * zswitch(ij) )
756               ! MV MP 2016
757               IF ( nn_pnd_scheme > 0 ) THEN
758                  v_ip (ji,jj,jl2) = v_ip (ji,jj,jl2)  + ( vprdg (ij) * rn_fpondrdg * fvol(ij)  +   &
759                                                       &   vprft (ij) * rn_fpondrft * zswitch(ij) )
760                  a_ip (ji,jj,jl2) = a_ip(ji,jj,jl2)   + ( aprdg2(ij) * rn_fpondrdg * farea +       & 
761                                                       &   aprft2(ij) * rn_fpondrft * zswitch(ji) )
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      ! SIMIP diagnostics
780      diag_dmi_dyn(:,:) = - wfx_dyn(:,:)     + rhoic * diag_trp_vi(:,:)
781      diag_dms_dyn(:,:) = - wfx_snw_dyn(:,:) + rhosn * diag_trp_vs(:,:)     
782      !
783   END SUBROUTINE ice_rdgrft_ridgeshift
784
785   SUBROUTINE ice_rdgrft_icestrength( kstrngth )
786      !!----------------------------------------------------------------------
787      !!                ***  ROUTINE ice_rdgrft_icestrength ***
788      !!
789      !! ** Purpose :   computes ice strength used in dynamics routines of ice thickness
790      !!
791      !! ** Method  :   Compute the strength of the ice pack, defined as the energy (J m-2)
792      !!              dissipated per unit area removed from the ice pack under compression,
793      !!              and assumed proportional to the change in potential energy caused
794      !!              by ridging. Note that only Hibler's formulation is stable and that
795      !!              ice strength has to be smoothed
796      !!
797      !! ** Inputs / Ouputs : kstrngth (what kind of ice strength we are using)
798      !!----------------------------------------------------------------------
799      INTEGER, INTENT(in) ::   kstrngth    ! = 1 for Rothrock formulation, 0 for Hibler (1979)
800      INTEGER             ::   ji,jj, jl   ! dummy loop indices
801      INTEGER             ::   ksmooth     ! smoothing the resistance to deformation
802      INTEGER             ::   numts_rm    ! number of time steps for the P smoothing
803      REAL(wp)            ::   zp, z1_3    ! local scalars
804      REAL(wp), DIMENSION(jpi,jpj) ::   zworka           ! temporary array used here
805      REAL(wp), DIMENSION(jpi,jpj) ::   zstrp1, zstrp2   ! strength at previous time steps
806      !!----------------------------------------------------------------------
807
808      !------------------------------------------------------------------------------!
809      ! 1) Initialize
810      !------------------------------------------------------------------------------!
811      strength(:,:) = 0._wp
812
813      !------------------------------------------------------------------------------!
814      ! 2) Compute thickness distribution of ridging and ridged ice
815      !------------------------------------------------------------------------------!
816      CALL ice_rdgrft_ridgeprep
817
818      !------------------------------------------------------------------------------!
819      ! 3) Rothrock(1975)'s method
820      !------------------------------------------------------------------------------!
821      IF( kstrngth == 1 ) THEN
822         z1_3 = 1._wp / 3._wp
823         DO jl = 1, jpl
824            DO jj= 1, jpj
825               DO ji = 1, jpi
826                  !
827                  IF( athorn(ji,jj,jl) > 0._wp ) THEN
828                     !----------------------------
829                     ! PE loss from deforming ice
830                     !----------------------------
831                     strength(ji,jj) = strength(ji,jj) - athorn(ji,jj,jl) * ht_i(ji,jj,jl) * ht_i(ji,jj,jl)
832
833                     !--------------------------
834                     ! PE gain from rafting ice
835                     !--------------------------
836                     strength(ji,jj) = strength(ji,jj) + 2._wp * araft(ji,jj,jl) * ht_i(ji,jj,jl) * ht_i(ji,jj,jl)
837
838                     !----------------------------
839                     ! PE gain from ridging ice
840                     !----------------------------
841                     strength(ji,jj) = strength(ji,jj) + aridge(ji,jj,jl) * krdg(ji,jj,jl) * z1_3 *  &
842                        &                              ( hrmax(ji,jj,jl) * hrmax(ji,jj,jl) +         &
843                        &                                hrmin(ji,jj,jl) * hrmin(ji,jj,jl) +         &
844                        &                                hrmax(ji,jj,jl) * hrmin(ji,jj,jl) ) 
845                        !!(a**3-b**3)/(a-b) = a*a+ab+b*b                     
846                  ENDIF
847                  !
848               END DO
849            END DO
850         END DO
851   
852         strength(:,:) = rn_pe_rdg * Cp * strength(:,:) / aksum(:,:) * tmask(:,:,1)
853                         ! where Cp = (g/2)*(rhow-rhoi)*(rhoi/rhow) and rn_pe_rdg accounts for frictional dissipation
854         ksmooth = 1
855
856      !------------------------------------------------------------------------------!
857      ! 4) Hibler (1979)' method
858      !------------------------------------------------------------------------------!
859      ELSE                      ! kstrngth ne 1:  Hibler (1979) form
860         !
861         strength(:,:) = rn_pstar * vt_i(:,:) * EXP( - rn_crhg * ( 1._wp - at_i(:,:) )  ) * tmask(:,:,1)
862         !
863         ksmooth = 1
864         !
865      ENDIF                     ! kstrngth
866      !
867      !------------------------------------------------------------------------------!
868      ! 5) Impact of brine volume
869      !------------------------------------------------------------------------------!
870      ! CAN BE REMOVED
871      IF( ln_icestr_bvf ) THEN
872         DO jj = 1, jpj
873            DO ji = 1, jpi
874               strength(ji,jj) = strength(ji,jj) * exp(-5.88*SQRT(MAX(bvm_i(ji,jj),0.0)))
875            END DO
876         END DO
877      ENDIF
878      !
879      !------------------------------------------------------------------------------!
880      ! 6) Smoothing ice strength
881      !------------------------------------------------------------------------------!
882      !
883      !-------------------
884      ! Spatial smoothing
885      !-------------------
886      IF ( ksmooth == 1 ) THEN
887
888         DO jj = 2, jpjm1
889            DO ji = 2, jpim1
890               IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp ) THEN
891                  zworka(ji,jj) = ( 4.0 * strength(ji,jj)              &
892                     &                  + strength(ji-1,jj) * tmask(ji-1,jj,1) + strength(ji+1,jj) * tmask(ji+1,jj,1) & 
893                     &                  + strength(ji,jj-1) * tmask(ji,jj-1,1) + strength(ji,jj+1) * tmask(ji,jj+1,1) &
894                     &            ) / ( 4.0 + tmask(ji-1,jj,1) + tmask(ji+1,jj,1) + tmask(ji,jj-1,1) + tmask(ji,jj+1,1) )
895               ELSE
896                  zworka(ji,jj) = 0._wp
897               ENDIF
898            END DO
899         END DO
900
901         DO jj = 2, jpjm1
902            DO ji = 2, jpim1
903               strength(ji,jj) = zworka(ji,jj)
904            END DO
905         END DO
906         CALL lbc_lnk( strength, 'T', 1. )
907
908      ENDIF
909
910      !--------------------
911      ! Temporal smoothing
912      !--------------------
913      IF ( ksmooth == 2 ) THEN
914
915         IF ( kt_ice == nit000 ) THEN
916            zstrp1(:,:) = 0._wp
917            zstrp2(:,:) = 0._wp
918         ENDIF
919
920         DO jj = 2, jpjm1
921            DO ji = 2, jpim1
922               IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp ) THEN
923                  numts_rm = 1 ! number of time steps for the running mean
924                  IF ( zstrp1(ji,jj) > 0._wp ) numts_rm = numts_rm + 1
925                  IF ( zstrp2(ji,jj) > 0._wp ) numts_rm = numts_rm + 1
926                  zp = ( strength(ji,jj) + zstrp1(ji,jj) + zstrp2(ji,jj) ) / numts_rm
927                  zstrp2(ji,jj) = zstrp1(ji,jj)
928                  zstrp1(ji,jj) = strength(ji,jj)
929                  strength(ji,jj) = zp
930               ENDIF
931            END DO
932         END DO
933
934         CALL lbc_lnk( strength, 'T', 1. )      ! Boundary conditions
935
936      ENDIF ! ksmooth
937      !
938   END SUBROUTINE ice_rdgrft_icestrength
939
940   SUBROUTINE ice_rdgrft_init
941      !!-------------------------------------------------------------------
942      !!                   ***  ROUTINE ice_rdgrft_init ***
943      !!
944      !! ** Purpose :   Physical constants and parameters linked
945      !!                to the mechanical ice redistribution
946      !!
947      !! ** Method  :   Read the namiceitdme namelist
948      !!                and check the parameters values
949      !!                called at the first timestep (nit000)
950      !!
951      !! ** input   :   Namelist namiceitdme
952      !!-------------------------------------------------------------------
953      INTEGER :: ios                 ! Local integer output status for namelist read
954      NAMELIST/namiceitdme/ rn_cs, nn_partfun, rn_gstar, rn_astar,             & 
955        &                   ln_ridging, rn_hstar, rn_por_rdg, rn_fsnowrdg, rn_fpondrdg, & 
956                            ln_rafting, rn_hraft, rn_craft,   rn_fsnowrft, rn_fpondrft
957      !!-------------------------------------------------------------------
958      !
959      REWIND( numnam_ice_ref )              ! Namelist namicetdme in reference namelist : Ice mechanical ice redistribution
960      READ  ( numnam_ice_ref, namiceitdme, IOSTAT = ios, ERR = 901)
961901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitdme in reference namelist', lwp )
962
963      REWIND( numnam_ice_cfg )              ! Namelist namiceitdme in configuration namelist : Ice mechanical ice redistribution
964      READ  ( numnam_ice_cfg, namiceitdme, IOSTAT = ios, ERR = 902 )
965902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitdme in configuration namelist', lwp )
966      IF(lwm) WRITE ( numoni, namiceitdme )
967      !
968      IF (lwp) THEN                          ! control print
969         WRITE(numout,*)
970         WRITE(numout,*)'ice_rdgrft_init : ice parameters for mechanical ice redistribution '
971         WRITE(numout,*)'~~~~~~~~~~~~~~~'
972         WRITE(numout,*)'   Fraction of shear energy contributing to ridging        rn_cs       = ', rn_cs 
973         WRITE(numout,*)'   Switch for part. function (0) linear (1) exponential    nn_partfun  = ', nn_partfun
974         WRITE(numout,*)'   Fraction of total ice coverage contributing to ridging  rn_gstar    = ', rn_gstar
975         WRITE(numout,*)'   Equivalent to G* for an exponential part function       rn_astar    = ', rn_astar
976         WRITE(numout,*)'   Ridging of ice sheets or not                            ln_ridging  = ', ln_ridging
977         WRITE(numout,*)'   Quantity playing a role in max ridged ice thickness     rn_hstar    = ', rn_hstar
978         WRITE(numout,*)'   Initial porosity of ridges                              rn_por_rdg  = ', rn_por_rdg
979         WRITE(numout,*)'   Fraction of snow volume conserved during ridging        rn_fsnowrdg = ', rn_fsnowrdg 
980         WRITE(numout,*)'   Fraction of pond volume conserved during ridging        rn_fpondrdg = ', rn_fpondrdg 
981         WRITE(numout,*)'   Rafting of ice sheets or not                            ln_rafting  = ', ln_rafting
982         WRITE(numout,*)'   Parmeter thickness (threshold between ridge-raft)       rn_hraft    = ', rn_hraft
983         WRITE(numout,*)'   Rafting hyperbolic tangent coefficient                  rn_craft    = ', rn_craft 
984         WRITE(numout,*)'   Fraction of snow volume conserved during ridging        rn_fsnowrft = ', rn_fsnowrft 
985         WRITE(numout,*)'   Fraction of pond volume conserved during rafting        rn_fpondrft = ', rn_fpondrft 
986      ENDIF
987      !
988   END SUBROUTINE ice_rdgrft_init
989
990#else
991   !!----------------------------------------------------------------------
992   !!   Default option         Empty module          NO LIM-3 sea-ice model
993   !!----------------------------------------------------------------------
994CONTAINS
995   SUBROUTINE ice_rdgrft           ! Empty routines
996   END SUBROUTINE ice_rdgrft
997   SUBROUTINE ice_rdgrft_icestrength
998   END SUBROUTINE ice_rdgrft_icestrength
999   SUBROUTINE ice_rdgrft_init
1000   END SUBROUTINE ice_rdgrft_init
1001#endif
1002   !!======================================================================
1003END MODULE icerdgrft
Note: See TracBrowser for help on using the repository browser.