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 @ 8486

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

changes in style - part1 - (now the code looks better txs to Gurvan's comments)

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