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

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

changes in style - part6 - commits of the day

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