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

Last change on this file since 8531 was 8531, checked in by clem, 3 years ago

changes in style - part6 - more clarity (still not finished)

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