source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 @ 7753

Last change on this file since 7753 was 7753, checked in by mocavero, 4 years ago

Reverting trunk to remove OpenMP

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