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

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

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

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

changes in style - part5 - continue changing init routines

File size: 48.7 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_icestrength   ! called by icerhg_evp
35   PUBLIC   ice_rdgrft_init          ! called by ice_stp
36   PUBLIC   ice_rdgrft_alloc         ! called by ice_init
37
38   ! Variables shared among ridging subroutines
39   REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   asum     ! sum of total ice and open water area
40   REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   aksum    ! ratio of area removed to area ridged
41   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   athorn   ! participation function; fraction of ridging/closing associated w/ category n
42   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   hrmin    ! minimum ridge thickness
43   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   hrmax    ! maximum ridge thickness
44   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   hraft    ! thickness of rafted ice
45   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   krdg     ! thickness of ridging ice / mean ridge thickness
46   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   aridge   ! participating ice ridging
47   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   araft    ! participating ice rafting
48   !
49   REAL(wp), PARAMETER ::   krdgmin = 1.1_wp    ! min ridge thickness multiplier
50   REAL(wp), PARAMETER ::   kraft   = 0.5_wp    ! rafting multipliyer
51   REAL(wp)            ::   zdrho               !
52   !
53   ! ** namelist (namiceitdme) **
54   REAL(wp) ::   rn_cs            !: fraction of shearing energy contributing to ridging           
55   INTEGER  ::   nn_partfun       !: participation function: =0 Thorndike et al. (1975), =1 Lipscomb et al. (2007)
56   REAL(wp) ::   rn_gstar         !: fractional area of young ice contributing to ridging
57   REAL(wp) ::   rn_astar         !: equivalent of G* for an exponential participation function
58   LOGICAL  ::   ln_ridging       !: ridging of ice or not                       
59   REAL(wp) ::   rn_hstar         !: thickness that determines the maximal thickness of ridged ice
60   REAL(wp) ::   rn_por_rdg       !: initial porosity of ridges (0.3 regular value)
61   REAL(wp) ::   rn_fsnowrdg      !: fractional snow loss to the ocean during ridging
62   REAL(wp) ::   rn_fpondrdg      !: fractional melt pond loss to the ocean during ridging
63   LOGICAL  ::   ln_rafting       !: rafting of ice or not                       
64   REAL(wp) ::   rn_hraft         !: threshold thickness (m) for rafting / ridging
65   REAL(wp) ::   rn_craft         !: coefficient for smoothness of the hyperbolic tangent in rafting
66   REAL(wp) ::   rn_fsnowrft      !: fractional snow loss to the ocean during ridging
67   REAL(wp) ::   rn_fpondrft      !: fractional snow loss to the ocean during rafting
68   !
69   !!----------------------------------------------------------------------
70   !! NEMO/ICE 4.0 , NEMO Consortium (2017)
71   !! $Id: icerdgrft.F90 8378 2017-07-26 13:55:59Z clem $
72   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
73   !!----------------------------------------------------------------------
74CONTAINS
75
76   INTEGER FUNCTION ice_rdgrft_alloc()
77      !!---------------------------------------------------------------------!
78      !!                ***  ROUTINE ice_rdgrft_alloc ***
79      !!---------------------------------------------------------------------!
80      ALLOCATE( asum (jpi,jpj)     , athorn(jpi,jpj,0:jpl) , aksum (jpi,jpj)     ,     &
81         &      hrmin(jpi,jpj,jpl) , hraft (jpi,jpj,jpl)   , aridge(jpi,jpj,jpl) ,     &
82         &      hrmax(jpi,jpj,jpl) , krdg  (jpi,jpj,jpl)   , araft (jpi,jpj,jpl) , STAT=ice_rdgrft_alloc )
83         !
84      IF( ice_rdgrft_alloc /= 0 )   CALL ctl_warn( 'ice_rdgrft_alloc: failed to allocate arrays' )
85      !
86   END FUNCTION ice_rdgrft_alloc
87
88
89   SUBROUTINE ice_rdgrft( kt )
90      !!---------------------------------------------------------------------!
91      !!                ***  ROUTINE ice_rdgrft ***
92      !!
93      !! ** Purpose :   computes the mechanical redistribution of ice thickness
94      !!
95      !! ** Method  :   Steps :
96      !!       1) Thickness categories boundaries, ice / o.w. concentrations
97      !!          Ridge preparation
98      !!       2) Dynamical inputs (closing rate, divu_adv, opning)
99      !!       3) Ridging iteration
100      !!       4) Ridging diagnostics
101      !!       5) Heat, salt and freshwater fluxes
102      !!       6) Compute increments of tate variables and come back to old values
103      !!
104      !! References :   Flato, G. M., and W. D. Hibler III, 1995, JGR, 100, 18,611-18,626.
105      !!                Hibler, W. D. III, 1980, MWR, 108, 1943-1973, 1980.
106      !!                Rothrock, D. A., 1975: JGR, 80, 4514-4519.
107      !!                Thorndike et al., 1975, JGR, 80, 4501-4513.
108      !!                Bitz et al., JGR, 2001
109      !!                Amundrud and Melling, JGR 2005
110      !!                Babko et al., JGR 2002
111      !!
112      !!     This routine is based on CICE code and authors William H. Lipscomb,
113      !!  and Elizabeth C. Hunke, LANL are gratefully acknowledged
114      !!--------------------------------------------------------------------!
115      INTEGER, INTENT(in) ::   kt     ! number of iteration
116      !!
117      INTEGER  ::   ji, jj, jk, jl     ! dummy loop index
118      INTEGER  ::   niter              ! local integer
119      INTEGER  ::   iterate_ridging    ! if =1, repeat the ridging
120      REAL(wp) ::   za                 ! local scalar
121      REAL(wp), DIMENSION(jpi,jpj) ::   closing_net     ! net rate at which area is removed    (1/s)
122      !                                                 ! (ridging ice area - area of new ridges) / dt
123      REAL(wp), DIMENSION(jpi,jpj) ::   divu_adv        ! divu as implied by transport scheme  (1/s)
124      REAL(wp), DIMENSION(jpi,jpj) ::   opning          ! rate of opening due to divergence/shear
125      REAL(wp), DIMENSION(jpi,jpj) ::   closing_gross   ! rate at which area removed, not counting area of new ridges
126      !
127      INTEGER, PARAMETER ::   nitermax = 20   
128      !
129      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
130      !!-----------------------------------------------------------------------------
131      IF( nn_timing == 1 )  CALL timing_start('icerdgrft')
132
133      IF( kt == nit000 ) THEN
134         IF(lwp) WRITE(numout,*)
135         IF(lwp) WRITE(numout,*)'icerdgrft : ice ridging and rafting'
136         IF(lwp) WRITE(numout,*)'~~~~~~~~~~'
137         !
138         zdrho = 0.5 * grav * (rau0-rhoic) * rhoic * r1_rau0      ! proport const for PE
139         !
140      ENDIF     
141      !                    ! conservation test
142      IF( ln_limdiachk )   CALL ice_cons_hsm(0, 'icerdgrft', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
143
144      !-----------------------------------------------------------------------------!
145      ! 1) Thickness categories boundaries, ice / o.w. concentrations, init_ons
146      !-----------------------------------------------------------------------------!
147      !
148      CALL ice_rdgrft_ridgeprep                             ! prepare ridging
149      !
150      DO jj = 1, jpj                                        ! Initialize arrays.
151         DO ji = 1, jpi
152            !-----------------------------------------------------------------------------!
153            ! 2) Dynamical inputs (closing rate, divu_adv, opning)
154            !-----------------------------------------------------------------------------!
155            !
156            ! 2.1 closing_net
157            !-----------------
158            ! Compute the net rate of closing due to convergence
159            ! and shear, based on Flato and Hibler (1995).
160            !
161            ! The energy dissipation rate is equal to the net closing rate
162            ! times the ice strength.
163            !
164            ! NOTE: The NET closing rate is equal to the rate that open water
165            !  area is removed, plus the rate at which ice area is removed by
166            !  ridging, minus the rate at which area is added in new ridges.
167            !  The GROSS closing rate is equal to the first two terms (open
168            !  water closing and thin ice ridging) without the third term
169            !  (thick, newly ridged ice).
170
171            closing_net(ji,jj) = rn_cs * 0.5_wp * ( delta_i(ji,jj) - ABS( divu_i(ji,jj) ) ) - MIN( divu_i(ji,jj), 0._wp )
172
173            ! 2.2 divu_adv
174            !--------------
175            ! Compute divu_adv, the divergence rate given by the transport/
176            ! advection scheme, which may not be equal to divu as computed
177            ! from the velocity field.
178            !
179            ! If divu_adv < 0, make sure the closing rate is large enough
180            ! to give asum = 1.0 after ridging.
181           
182            divu_adv(ji,jj) = ( 1._wp - asum(ji,jj) ) * r1_rdtice  ! asum found in ridgeprep
183
184            IF( divu_adv(ji,jj) < 0._wp )   closing_net(ji,jj) = MAX( closing_net(ji,jj), -divu_adv(ji,jj) )
185
186            ! 2.3 opning
187            !------------
188            ! Compute the (non-negative) opening rate that will give asum = 1.0 after ridging.
189            opning(ji,jj) = closing_net(ji,jj) + divu_adv(ji,jj)
190         END DO
191      END DO
192
193      !-----------------------------------------------------------------------------!
194      ! 3) Ridging iteration
195      !-----------------------------------------------------------------------------!
196      niter           = 1                 ! iteration counter
197      iterate_ridging = 1
198
199      DO WHILE ( iterate_ridging > 0 .AND. niter < nitermax )
200
201         ! 3.2 closing_gross
202         !-----------------------------------------------------------------------------!
203         ! Based on the ITD of ridging and ridged ice, convert the net
204         !  closing rate to a gross closing rate. 
205         ! NOTE: 0 < aksum <= 1
206         closing_gross(:,:) = closing_net(:,:) / aksum(:,:)
207
208         ! correction to closing rate and opening if closing rate is excessive
209         !---------------------------------------------------------------------
210         ! Reduce the closing rate if more than 100% of the open water
211         ! would be removed.  Reduce the opening rate proportionately.
212         DO jj = 1, jpj
213            DO ji = 1, jpi
214               za   = ( opning(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) ) * rdt_ice
215               IF    ( za < 0._wp .AND. za > - ato_i(ji,jj) ) THEN                  ! would lead to negative ato_i
216                  opning(ji,jj) = athorn(ji,jj,0) * closing_gross(ji,jj) - ato_i(ji,jj) * r1_rdtice 
217               ELSEIF( za > 0._wp .AND. za > ( asum(ji,jj) - ato_i(ji,jj) ) ) THEN  ! would lead to ato_i > asum
218                  opning(ji,jj) = athorn(ji,jj,0) * closing_gross(ji,jj) + ( asum(ji,jj) - ato_i(ji,jj) ) * r1_rdtice 
219               ENDIF
220            END DO
221         END DO
222
223         ! correction to closing rate / opening if excessive ice removal
224         !---------------------------------------------------------------
225         ! Reduce the closing rate if more than 100% of any ice category
226         ! would be removed.  Reduce the opening rate proportionately.
227         DO jl = 1, jpl
228            DO jj = 1, jpj
229               DO ji = 1, jpi
230                  za = athorn(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice
231                  IF( za  >  a_i(ji,jj,jl) ) THEN
232                     closing_gross(ji,jj) = closing_gross(ji,jj) * a_i(ji,jj,jl) / za
233                  ENDIF
234               END DO
235            END DO
236         END DO
237
238         ! 3.3 Redistribute area, volume, and energy.
239         !-----------------------------------------------------------------------------!
240         CALL ice_rdgrft_ridgeshift( opning, closing_gross )
241         
242         ! 3.4 Compute total area of ice plus open water after ridging.
243         !-----------------------------------------------------------------------------!
244         ! This is in general not equal to one because of divergence during transport
245         asum(:,:) = ato_i(:,:) + SUM( a_i, dim=3 )
246
247         ! 3.5 Do we keep on iterating?
248         !-----------------------------------------------------------------------------!
249         ! Check whether asum = 1.  If not (because the closing and opening
250         ! rates were reduced above), ridge again with new rates.
251         iterate_ridging = 0
252         DO jj = 1, jpj
253            DO ji = 1, jpi
254               IF( ABS( asum(ji,jj) - 1._wp ) < epsi10 ) THEN
255                  closing_net(ji,jj) = 0._wp
256                  opning     (ji,jj) = 0._wp
257                  ato_i      (ji,jj) = MAX( 0._wp, 1._wp - SUM( a_i(ji,jj,:) ) )
258               ELSE
259                  iterate_ridging    = 1
260                  divu_adv   (ji,jj) = ( 1._wp - asum(ji,jj) ) * r1_rdtice
261                  closing_net(ji,jj) = MAX( 0._wp, -divu_adv(ji,jj) )
262                  opning     (ji,jj) = MAX( 0._wp,  divu_adv(ji,jj) )
263               ENDIF
264            END DO
265         END DO
266         IF( lk_mpp )   CALL mpp_max( iterate_ridging )
267
268         ! Repeat if necessary.
269         ! NOTE: If strength smoothing is turned on, the ridging must be
270         ! iterated globally because of the boundary update in the smoothing.
271         niter = niter + 1
272         !
273         IF( iterate_ridging == 1 ) THEN
274            CALL ice_rdgrft_ridgeprep
275            IF( niter  >  nitermax ) THEN
276               WRITE(numout,*) ' ALERTE : non-converging ridging scheme '
277               WRITE(numout,*) ' niter, iterate_ridging ', niter, iterate_ridging
278            ENDIF
279         ENDIF
280         !
281      END DO !! on the do while over iter
282
283      CALL ice_var_agg( 1 ) 
284
285      !-----------------------------------------------------------------------------!
286      ! control prints
287      !-----------------------------------------------------------------------------!
288      !                    ! conservation test
289      IF( ln_limdiachk )   CALL ice_cons_hsm(1, 'icerdgrft', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
290
291      !                    ! control prints
292      IF( ln_ctl       )   CALL ice_prt3D( 'icerdgrft' )
293      !
294      IF( nn_timing == 1 )  CALL timing_stop('icerdgrft')
295      !
296   END SUBROUTINE ice_rdgrft
297
298
299   SUBROUTINE ice_rdgrft_ridgeprep
300      !!---------------------------------------------------------------------!
301      !!                ***  ROUTINE ice_rdgrft_ridgeprep ***
302      !!
303      !! ** Purpose :   preparation for ridging and strength calculations
304      !!
305      !! ** Method  :   Compute the thickness distribution of the ice and open water
306      !!              participating in ridging and of the resulting ridges.
307      !!---------------------------------------------------------------------!
308      INTEGER  ::   ji, jj, jl    ! dummy loop indices
309      REAL(wp) ::   z1_gstar, z1_astar, zhmean, zdummy   ! local scalar
310      REAL(wp), DIMENSION(jpi,jpj,-1:jpl) ::   zGsum     ! zGsum(n) = sum of areas in categories 0 to n
311      !------------------------------------------------------------------------------!
312
313      z1_gstar      = 1._wp / rn_gstar   
314      z1_astar      = 1._wp / rn_astar   
315
316      CALL ice_var_zapsmall   ! Zero out categories with very small areas
317
318      !                       ! Ice thickness needed for rafting
319      WHERE( a_i(:,:,:) >= epsi20 )   ;   ht_i(:,:,:) = v_i (:,:,:) / a_i(:,:,:)
320      ELSEWHERE                       ;   ht_i(:,:,:) = 0._wp
321      END WHERE
322
323      !------------------------------------------------------------------------------!
324      ! 1) Participation function
325      !------------------------------------------------------------------------------!
326      !
327      ! Compute total area of ice plus open water.
328      ! This is in general not equal to one because of divergence during transport
329      asum(:,:) = ato_i(:,:) + SUM( a_i, dim=3 )
330      !
331      ! Compute cumulative thickness distribution function
332      ! Compute the cumulative thickness distribution function zGsum,
333      ! where zGsum(n) is the fractional area in categories 0 to n.
334      ! initial value (in h = 0) equals open water area
335      zGsum(:,:,-1) = 0._wp
336      zGsum(:,:,0 ) = ato_i(:,:) / asum(:,:)
337      DO jl = 1, jpl
338         zGsum(:,:,jl) = ( ato_i(:,:) + SUM( a_i(:,:,1:jl), dim=3 ) ) / asum(:,:) ! sum(1:jl) is ok (and not jpl)
339      END DO
340
341      ! 1.3 Compute participation function a(h) = b(h).g(h) (athorn)
342      !--------------------------------------------------------------------------------------------------
343      ! Compute the participation function athorn; this is analogous to
344      ! a(h) = b(h)g(h) as defined in Thorndike et al. (1975).
345      ! area lost from category n due to ridging/closing
346      ! athorn(n)   = total area lost due to ridging/closing
347      ! assume b(h) = (2/Gstar) * (1 - G(h)/Gstar).
348      !
349      ! The expressions for athorn are found by integrating b(h)g(h) between
350      ! the category boundaries.
351      ! athorn is always >= 0 and SUM(athorn(0:jpl))=1
352      !-----------------------------------------------------------------
353      !
354      IF( nn_partfun == 0 ) THEN       !--- Linear formulation (Thorndike et al., 1975)
355         DO jl = 0, jpl   
356            DO jj = 1, jpj 
357               DO ji = 1, jpi
358                  IF    ( zGsum(ji,jj,jl)   < rn_gstar ) THEN
359                     athorn(ji,jj,jl) = z1_gstar * ( zGsum(ji,jj,jl) - zGsum(ji,jj,jl-1) ) * &
360                        &                          ( 2._wp - ( zGsum(ji,jj,jl-1) + zGsum(ji,jj,jl) ) * z1_gstar )
361                  ELSEIF( zGsum(ji,jj,jl-1) < rn_gstar ) THEN
362                     athorn(ji,jj,jl) = z1_gstar * ( rn_gstar        - zGsum(ji,jj,jl-1) ) *  &
363                        &                          ( 2._wp - ( zGsum(ji,jj,jl-1) + rn_gstar        ) * z1_gstar )
364                  ELSE
365                     athorn(ji,jj,jl) = 0._wp
366                  ENDIF
367               END DO
368            END DO
369         END DO
370         !
371      ELSE                             !--- Exponential, more stable formulation (Lipscomb et al, 2007)
372         !                       
373         zdummy = 1._wp / ( 1._wp - EXP(-z1_astar) )        ! precompute exponential terms using zGsum as a work array
374         DO jl = -1, jpl
375            zGsum(:,:,jl) = EXP( -zGsum(:,:,jl) * z1_astar ) * zdummy
376         END DO
377         DO jl = 0, jpl
378            athorn(:,:,jl) = zGsum(:,:,jl-1) - zGsum(:,:,jl)
379         END DO
380         !
381      ENDIF
382
383      !                                !--- Ridging and rafting participation concentrations
384      IF( ln_rafting .AND. ln_ridging ) THEN             !- ridging & rafting
385         DO jl = 1, jpl
386            DO jj = 1, jpj 
387               DO ji = 1, jpi
388                  aridge(ji,jj,jl) = ( 1._wp + TANH ( rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) ) ) * 0.5_wp * athorn(ji,jj,jl)
389                  araft (ji,jj,jl) = athorn(ji,jj,jl) - aridge(ji,jj,jl)
390               END DO
391            END DO
392         END DO
393      ELSEIF( ln_ridging .AND. .NOT. ln_rafting ) THEN   !- ridging alone
394         aridge(:,:,:) = athorn(:,:,1:jpl)
395         araft (:,:,:) = 0._wp
396      ELSEIF( ln_rafting .AND. .NOT. ln_ridging ) THEN   !- rafting alone   
397        aridge(:,:,:) = 0._wp
398         araft (:,:,:) = athorn(:,:,1:jpl)
399      ELSE                                               !- no ridging & no rafting
400         aridge(:,:,:) = 0._wp
401         araft (:,:,:) = 0._wp         
402      ENDIF
403
404      !-----------------------------------------------------------------
405      ! 2) Transfer function
406      !-----------------------------------------------------------------
407      ! Compute max and min ridged ice thickness for each ridging category.
408      ! Assume ridged ice is uniformly distributed between hrmin and hrmax.
409      !
410      ! This parameterization is a modified version of Hibler (1980).
411      ! The mean ridging thickness, zhmean, is proportional to hi^(0.5)
412      !  and for very thick ridging ice must be >= krdgmin*hi
413      !
414      ! The minimum ridging thickness, hrmin, is equal to 2*hi
415      !  (i.e., rafting) and for very thick ridging ice is
416      !  constrained by hrmin <= (zhmean + hi)/2.
417      !
418      ! The maximum ridging thickness, hrmax, is determined by
419      !  zhmean and hrmin.
420      !
421      ! These modifications have the effect of reducing the ice strength
422      ! (relative to the Hibler formulation) when very thick ice is
423      ! ridging.
424      !
425      ! aksum = net area removed/ total area removed
426      ! where total area removed = area of ice that ridges
427      !         net area removed = total area removed - area of new ridges
428      !-----------------------------------------------------------------
429
430      aksum(:,:) = athorn(:,:,0)
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 ( athorn(ji,jj,jl) > 0._wp ) THEN
436                  zhmean          = MAX( SQRT( rn_hstar * ht_i(ji,jj,jl) ), ht_i(ji,jj,jl) * krdgmin )
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) / kraft
440                  krdg (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 - krdg(ji,jj,jl) )    &
444                     &                        + araft (ji,jj,jl) * ( 1._wp - kraft          )
445               ELSE
446                  hrmin(ji,jj,jl)  = 0._wp 
447                  hrmax(ji,jj,jl)  = 0._wp 
448                  hraft(ji,jj,jl)  = 0._wp 
449                  krdg (ji,jj,jl)  = 1._wp
450               ENDIF
451            END DO
452         END DO
453      END DO
454      !
455   END SUBROUTINE ice_rdgrft_ridgeprep
456
457
458   SUBROUTINE ice_rdgrft_ridgeshift( opning, closing_gross )
459      !!----------------------------------------------------------------------
460      !!                ***  ROUTINE ice_rdgrft_icestrength ***
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) - athorn(ji,jj,0) * closing_gross(ji,jj) ) * rdt_ice )
520         END DO
521      END DO
522
523      !-----------------------------------------------------------------
524      ! 3) 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         ! 3.1) Identify grid cells with nonzero ridging
533         !------------------------------------------------
534         icells = 0
535         DO jj = 1, jpj
536            DO ji = 1, jpi
537               IF( athorn(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            ! 3.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            ! 3.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) * krdg(ji,jj,jl1)
560            arft2(ij) = arft1(ij) * kraft
561
562            !--------------------------------------------------------------------------
563            ! 3.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_por_rdg )
568            vsw  (ij) = vrdg1(ij) * rn_por_rdg
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) * krdg(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) * krdg(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) * kraft
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) * kraft 
598
599            !-----------------------------------------------------------------
600            ! 3.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            ! 3.7 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_fsnowrdg )   & 
624               &                                      + rhosn * vsrft(ij) * ( 1._wp - rn_fsnowrft ) ) * r1_rdtice  ! fresh water source for ocean
625
626            hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ( - esrdg(ij) * ( 1._wp - rn_fsnowrdg )         & 
627               &                                - esrft(ij) * ( 1._wp - rn_fsnowrft ) ) * r1_rdtice        ! heat sink for ocean (<0, W.m-2)
628
629            ! MV MP 2016
630            !------------------------------------------           
631            ! 3.X 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_fpondrdg )   & 
636               &                                 + rhofw * vprft(ij) * ( 1._wp - rn_fpondrft ) ) * r1_rdtice  ! fresh water source for ocean
637            ENDIF
638            ! END MV MP 2016
639
640            !-----------------------------------------------------------------
641            ! 3.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         ! 3.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         ! 4) 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_fsnowrdg * fvol(ij)  +  &
722                  &                                        vsrft (ij) * rn_fsnowrft * zswitch(ij) )
723               e_s  (ji,jj,1,jl2) = e_s  (ji,jj,1,jl2) + ( esrdg (ij) * rn_fsnowrdg * fvol(ij)  +  &
724                  &                                        esrft (ij) * rn_fsnowrft * 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_fpondrdg * fvol   (ij)   &
728                     &                                   + vprft (ij) * rn_fpondrft * zswitch(ij)   )
729                  a_ip (ji,jj,jl2) = a_ip(ji,jj,jl2) + (   aprdg2(ij) * rn_fpondrdg * farea         & 
730                     &                                   + aprft2(ij) * rn_fpondrft * 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 ice_rdgrft_ridgeshift
748
749
750   SUBROUTINE ice_rdgrft_icestrength( kstrngth )
751      !!----------------------------------------------------------------------
752      !!                ***  ROUTINE ice_rdgrft_icestrength ***
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      !! ** Inputs / Ouputs : kstrngth (what kind of ice strength we are using)
763      !!----------------------------------------------------------------------
764      INTEGER, INTENT(in) ::   kstrngth    ! = 1 for Rothrock formulation, 0 for Hibler (1979)
765      !
766      INTEGER             ::   ji,jj, jl   ! dummy loop indices
767      INTEGER             ::   ismooth     ! smoothing the resistance to deformation
768      INTEGER             ::   itframe     ! number of time steps for the P smoothing
769      REAL(wp)            ::   zp, z1_3    ! local scalars
770      REAL(wp), DIMENSION(jpi,jpj) ::   zworka           ! temporary array used here
771      REAL(wp), DIMENSION(jpi,jpj) ::   zstrp1, zstrp2   ! strength at previous time steps
772      !!----------------------------------------------------------------------
773
774      !                              !--------------------------------------------------!
775      CALL ice_rdgrft_ridgeprep      ! Thickness distribution of ridging and ridged ice !
776      !                              !--------------------------------------------------!
777
778      !                              !--------------------------------------------------!
779      IF( kstrngth == 1 ) THEN       ! Ice strength => Rothrock (1975) method           !
780      !                              !--------------------------------------------------!
781         z1_3 = 1._wp / 3._wp
782         DO jl = 1, jpl
783            WHERE( athorn(:,:,jl) > 0._wp )
784               strength(:,:) =  -         athorn(:,:,jl) * ht_i(:,:,jl) * ht_i(:,:,jl)   &  ! PE loss from deforming ice
785                  &             + 2._wp * araft (:,:,jl) * ht_i(:,:,jl) * ht_i(:,:,jl)   &  ! PE gain from rafting ice
786                  &             +         aridge(:,:,jl) * krdg(:,:,jl) * z1_3 *   &        ! PE gain from ridging ice
787                  &                      ( hrmax(:,:,jl) * hrmax(:,:,jl) +         &
788                  &                        hrmin(:,:,jl) * hrmin(:,:,jl) +         &
789                  &                        hrmax(:,:,jl) * hrmin(:,:,jl) )
790            ELSEWHERE
791               strength(:,:) = 0._wp
792            END WHERE
793         END DO
794         strength(:,:) = rn_pe_rdg * zdrho * strength(:,:) / aksum(:,:) * tmask(:,:,1)
795                         ! where zdrho = (g/2)*(rhow-rhoi)*(rhoi/rhow) and rn_pe_rdg accounts for frictional dissipation
796         ismooth = 1
797      !                              !--------------------------------------------------!
798      ELSE                           ! Ice strength => Hibler (1979) method             !
799      !                              !--------------------------------------------------!
800         strength(:,:) = rn_pstar * vt_i(:,:) * EXP( - rn_crhg * ( 1._wp - at_i(:,:) )  ) * tmask(:,:,1)
801         !
802         ismooth = 1
803         !
804      ENDIF
805      !                              !--------------------------------------------------!
806      IF( ln_icestr_bvf ) THEN       ! Impact of brine volume                           !
807      !                              !--------------------------------------------------!
808      ! CAN BE REMOVED
809         DO jj = 1, jpj
810            DO ji = 1, jpi
811               strength(ji,jj) = strength(ji,jj) * exp(-5.88*SQRT(MAX(bvm_i(ji,jj),0.0)))
812            END DO
813         END DO
814      ENDIF
815      !                              !--------------------------------------------------!
816      SELECT CASE( ismooth )         ! Smoothing ice strength                           !
817      !                              !--------------------------------------------------!
818      !
819      CASE( 1 )               !--- Spatial smoothing
820         DO jj = 2, jpjm1
821            DO ji = 2, jpim1
822               IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp ) THEN
823                  zworka(ji,jj) = ( 4.0 * strength(ji,jj)              &
824                     &                  + strength(ji-1,jj) * tmask(ji-1,jj,1) + strength(ji+1,jj) * tmask(ji+1,jj,1) & 
825                     &                  + strength(ji,jj-1) * tmask(ji,jj-1,1) + strength(ji,jj+1) * tmask(ji,jj+1,1) &
826                     &            ) / ( 4.0 + tmask(ji-1,jj,1) + tmask(ji+1,jj,1) + tmask(ji,jj-1,1) + tmask(ji,jj+1,1) )
827               ELSE
828                  zworka(ji,jj) = 0._wp
829               ENDIF
830            END DO
831         END DO
832         
833         DO jj = 2, jpjm1
834            DO ji = 2, jpim1
835               strength(ji,jj) = zworka(ji,jj)
836            END DO
837         END DO
838         CALL lbc_lnk( strength, 'T', 1. )
839         !
840      CASE( 2 )               !--- Temporal smoothing
841         IF ( kt_ice == nit000 ) THEN
842            zstrp1(:,:) = 0._wp
843            zstrp2(:,:) = 0._wp
844         ENDIF
845         !
846         DO jj = 2, jpjm1
847            DO ji = 2, jpim1
848               IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp ) THEN
849                  itframe = 1 ! number of time steps for the running mean
850                  IF ( zstrp1(ji,jj) > 0._wp ) itframe = itframe + 1
851                  IF ( zstrp2(ji,jj) > 0._wp ) itframe = itframe + 1
852                  zp = ( strength(ji,jj) + zstrp1(ji,jj) + zstrp2(ji,jj) ) / itframe
853                  zstrp2  (ji,jj) = zstrp1  (ji,jj)
854                  zstrp1  (ji,jj) = strength(ji,jj)
855                  strength(ji,jj) = zp
856               ENDIF
857            END DO
858         END DO
859         CALL lbc_lnk( strength, 'T', 1. )
860         !
861      END SELECT
862      !
863   END SUBROUTINE ice_rdgrft_icestrength
864
865
866   SUBROUTINE ice_rdgrft_init
867      !!-------------------------------------------------------------------
868      !!                  ***  ROUTINE ice_rdgrft_init ***
869      !!
870      !! ** Purpose :   Physical constants and parameters linked
871      !!                to the mechanical ice redistribution
872      !!
873      !! ** Method  :   Read the namiceitdme namelist
874      !!                and check the parameters values
875      !!                called at the first timestep (nit000)
876      !!
877      !! ** input   :   Namelist namiceitdme
878      !!-------------------------------------------------------------------
879      INTEGER :: ios                 ! Local integer output status for namelist read
880      !!
881      NAMELIST/namiceitdme/ rn_cs     , nn_partfun, rn_gstar  , rn_astar   ,                & 
882         &                  ln_ridging, rn_hstar  , rn_por_rdg, rn_fsnowrdg, rn_fpondrdg,   & 
883         &                  ln_rafting, rn_hraft  , rn_craft  , rn_fsnowrft, rn_fpondrft
884      !!-------------------------------------------------------------------
885      !
886      REWIND( numnam_ice_ref )              ! Namelist namicetdme in reference namelist : Ice mechanical ice redistribution
887      READ  ( numnam_ice_ref, namiceitdme, IOSTAT = ios, ERR = 901)
888901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitdme in reference namelist', lwp )
889      !
890      REWIND( numnam_ice_cfg )              ! Namelist namiceitdme in configuration namelist : Ice mechanical ice redistribution
891      READ  ( numnam_ice_cfg, namiceitdme, IOSTAT = ios, ERR = 902 )
892902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitdme in configuration namelist', lwp )
893      IF(lwm) WRITE ( numoni, namiceitdme )
894      !
895      IF (lwp) THEN                          ! control print
896         WRITE(numout,*)
897         WRITE(numout,*) 'ice_rdgrft_init : ice parameters for mechanical ice redistribution '
898         WRITE(numout,*) '~~~~~~~~~~~~~~~'
899         WRITE(numout,*) '   Namelist namiceitdme'
900         WRITE(numout,*) '      Fraction of shear energy contributing to ridging        rn_cs       = ', rn_cs 
901         WRITE(numout,*) '      Switch for part. function (0) linear (1) exponential    nn_partfun  = ', nn_partfun
902         WRITE(numout,*) '      Fraction of total ice coverage contributing to ridging  rn_gstar    = ', rn_gstar
903         WRITE(numout,*) '      Equivalent to G* for an exponential part function       rn_astar    = ', rn_astar
904         WRITE(numout,*) '      Ridging of ice sheets or not                            ln_ridging  = ', ln_ridging
905         WRITE(numout,*) '      Quantity playing a role in max ridged ice thickness     rn_hstar    = ', rn_hstar
906         WRITE(numout,*) '      Initial porosity of ridges                              rn_por_rdg  = ', rn_por_rdg
907         WRITE(numout,*) '      Fraction of snow volume conserved during ridging        rn_fsnowrdg = ', rn_fsnowrdg 
908         WRITE(numout,*) '      Fraction of pond volume conserved during ridging        rn_fpondrdg = ', rn_fpondrdg 
909         WRITE(numout,*) '      Rafting of ice sheets or not                            ln_rafting  = ', ln_rafting
910         WRITE(numout,*) '      Parmeter thickness (threshold between ridge-raft)       rn_hraft    = ', rn_hraft
911         WRITE(numout,*) '      Rafting hyperbolic tangent coefficient                  rn_craft    = ', rn_craft 
912         WRITE(numout,*) '      Fraction of snow volume conserved during ridging        rn_fsnowrft = ', rn_fsnowrft 
913         WRITE(numout,*) '      Fraction of pond volume conserved during rafting        rn_fpondrft = ', rn_fpondrft 
914      ENDIF
915      !
916   END SUBROUTINE ice_rdgrft_init
917
918#else
919   !!----------------------------------------------------------------------
920   !!   Default option         Empty module          NO LIM-3 sea-ice model
921   !!----------------------------------------------------------------------
922#endif
923
924   !!======================================================================
925END MODULE icerdgrft
Note: See TracBrowser for help on using the repository browser.