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

source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rdgrft.F90 @ 8534

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

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