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

Last change on this file since 5134 was 5134, checked in by clem, 6 years ago

LIM3: changes to constrain ice thickness after advection. Ice volume and concentration are advected while ice thickness is deduced. This could result in overly high thicknesses associated with very low concentrations (in the 5th category)

  • Property svn:keywords set to Id
File size: 61.7 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 dom_ice          ! LIM domain
21   USE limvar           ! LIM
22   USE lbclnk           ! lateral boundary condition - MPP exchanges
23   USE lib_mpp          ! MPP library
24   USE wrk_nemo         ! work arrays
25   USE prtctl           ! Print control
26
27   USE in_out_manager   ! I/O manager
28   USE iom              ! I/O manager
29   USE lib_fortran      ! glob_sum
30   USE limdiahsb
31   USE timing           ! Timing
32   USE limcons          ! conservation tests
33
34   IMPLICIT NONE
35   PRIVATE
36
37   PUBLIC   lim_itd_me               ! called by ice_stp
38   PUBLIC   lim_itd_me_icestrength
39   PUBLIC   lim_itd_me_init
40   PUBLIC   lim_itd_me_alloc        ! called by sbc_lim_init
41
42   !-----------------------------------------------------------------------
43   ! Variables shared among ridging subroutines
44   !-----------------------------------------------------------------------
45   REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   asum     ! sum of total ice and open water area
46   REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   aksum    ! ratio of area removed to area ridged
47   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   athorn   ! participation function; fraction of ridging/
48   !                                                     ! closing associated w/ category n
49   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   hrmin    ! minimum ridge thickness
50   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   hrmax    ! maximum ridge thickness
51   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   hraft    ! thickness of rafted ice
52   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   krdg     ! mean ridge thickness/thickness of ridging ice
53   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   aridge   ! participating ice ridging
54   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   araft    ! participating ice rafting
55
56   REAL(wp), PARAMETER ::   krdgmin = 1.1_wp    ! min ridge thickness multiplier
57   REAL(wp), PARAMETER ::   kraft   = 2.0_wp    ! rafting multipliyer
58   REAL(wp), PARAMETER ::   kamax   = 1.0_wp    ! max of ice area authorized (clem: scheme is not stable if kamax <= 0.99)
59
60   REAL(wp) ::   Cp                             !
61   !
62   !-----------------------------------------------------------------------
63   ! Ridging diagnostic arrays for history files
64   !-----------------------------------------------------------------------
65   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   dardg1dt   ! rate of fractional area loss by ridging ice (1/s)
66   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   dardg2dt   ! rate of fractional area gain by new ridges (1/s)
67   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   dvirdgdt   ! rate of ice volume ridged (m/s)
68   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   opening    ! rate of opening due to divergence/shear (1/s)
69   !
70   !!----------------------------------------------------------------------
71   !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010)
72   !! $Id$
73   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
74   !!----------------------------------------------------------------------
75CONTAINS
76
77   INTEGER FUNCTION lim_itd_me_alloc()
78      !!---------------------------------------------------------------------!
79      !!                ***  ROUTINE lim_itd_me_alloc ***
80      !!---------------------------------------------------------------------!
81      ALLOCATE(                                                                     &
82         !* Variables shared among ridging subroutines
83         &      asum (jpi,jpj)     , athorn(jpi,jpj,0:jpl)                    ,     &
84         &      aksum(jpi,jpj)                                                ,     &
85         !
86         &      hrmin(jpi,jpj,jpl) , hraft(jpi,jpj,jpl) , aridge(jpi,jpj,jpl) ,     &
87         &      hrmax(jpi,jpj,jpl) , krdg (jpi,jpj,jpl) , araft (jpi,jpj,jpl) ,     &
88         !
89         !* Ridging diagnostic arrays for history files
90         &      dardg1dt(jpi,jpj)  , dardg2dt(jpi,jpj)                        ,     & 
91         &      dvirdgdt(jpi,jpj)  , opening(jpi,jpj)                         , STAT=lim_itd_me_alloc )
92         !
93      IF( lim_itd_me_alloc /= 0 )   CALL ctl_warn( 'lim_itd_me_alloc: failed to allocate arrays' )
94      !
95   END FUNCTION lim_itd_me_alloc
96
97
98   SUBROUTINE lim_itd_me
99      !!---------------------------------------------------------------------!
100      !!                ***  ROUTINE lim_itd_me ***
101      !!
102      !! ** Purpose :   computes the mechanical redistribution of ice thickness
103      !!
104      !! ** Method  :   Steps :
105      !!       1) Thickness categories boundaries, ice / o.w. concentrations
106      !!          Ridge preparation
107      !!       2) Dynamical inputs (closing rate, divu_adv, opning)
108      !!       3) Ridging iteration
109      !!       4) Ridging diagnostics
110      !!       5) Heat, salt and freshwater fluxes
111      !!       6) Compute increments of tate variables and come back to old values
112      !!
113      !! References :   Flato, G. M., and W. D. Hibler III, 1995, JGR, 100, 18,611-18,626.
114      !!                Hibler, W. D. III, 1980, MWR, 108, 1943-1973, 1980.
115      !!                Rothrock, D. A., 1975: JGR, 80, 4514-4519.
116      !!                Thorndike et al., 1975, JGR, 80, 4501-4513.
117      !!                Bitz et al., JGR, 2001
118      !!                Amundrud and Melling, JGR 2005
119      !!                Babko et al., JGR 2002
120      !!
121      !!     This routine is based on CICE code and authors William H. Lipscomb,
122      !!  and Elizabeth C. Hunke, LANL are gratefully acknowledged
123      !!--------------------------------------------------------------------!
124      INTEGER  ::   ji, jj, jk, jl        ! dummy loop index
125      INTEGER  ::   niter                 ! local integer
126      INTEGER  ::   iterate_ridging       ! if true, repeat the ridging
127      REAL(wp) ::   za, zfac              ! local scalar
128      CHARACTER (len = 15) ::   fieldid
129      REAL(wp), POINTER, DIMENSION(:,:) ::   closing_net     ! net rate at which area is removed    (1/s)
130                                                             ! (ridging ice area - area of new ridges) / dt
131      REAL(wp), POINTER, DIMENSION(:,:) ::   divu_adv        ! divu as implied by transport scheme  (1/s)
132      REAL(wp), POINTER, DIMENSION(:,:) ::   opning          ! rate of opening due to divergence/shear
133      REAL(wp), POINTER, DIMENSION(:,:) ::   closing_gross   ! rate at which area removed, not counting area of new ridges
134      REAL(wp), POINTER, DIMENSION(:,:) ::   msnow_mlt       ! mass of snow added to ocean (kg m-2)
135      REAL(wp), POINTER, DIMENSION(:,:) ::   esnow_mlt       ! energy needed to melt snow in ocean (J m-2)
136      REAL(wp), POINTER, DIMENSION(:,:) ::   vt_i_init, vt_i_final  !  ice volume summed over categories
137      !
138      INTEGER, PARAMETER ::   nitermax = 20   
139      !
140      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
141      !!-----------------------------------------------------------------------------
142      IF( nn_timing == 1 )  CALL timing_start('limitd_me')
143
144      CALL wrk_alloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final )
145
146      IF(ln_ctl) THEN
147         CALL prt_ctl(tab2d_1=ato_i , clinfo1=' lim_itd_me: ato_i  : ', tab2d_2=at_i   , clinfo2=' at_i    : ')
148         CALL prt_ctl(tab2d_1=divu_i, clinfo1=' lim_itd_me: divu_i : ', tab2d_2=delta_i, clinfo2=' delta_i : ')
149      ENDIF
150
151      IF( ln_limdyn ) THEN          !   Start ridging and rafting   !
152
153      ! conservation test
154      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limitd_me', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
155
156      !-----------------------------------------------------------------------------!
157      ! 1) Thickness categories boundaries, ice / o.w. concentrations, init_ons
158      !-----------------------------------------------------------------------------!
159      Cp = 0.5 * grav * (rau0-rhoic) * rhoic * r1_rau0             ! proport const for PE
160      !
161      CALL lim_itd_me_ridgeprep                                    ! prepare ridging
162      !
163      IF( con_i)   CALL lim_column_sum( jpl, v_i, vt_i_init )      ! conservation check
164
165      DO jj = 1, jpj                                     ! Initialize arrays.
166         DO ji = 1, jpi
167            msnow_mlt(ji,jj) = 0._wp
168            esnow_mlt(ji,jj) = 0._wp
169            dardg1dt (ji,jj) = 0._wp
170            dardg2dt (ji,jj) = 0._wp
171            dvirdgdt (ji,jj) = 0._wp
172            opening  (ji,jj) = 0._wp
173
174            !-----------------------------------------------------------------------------!
175            ! 2) Dynamical inputs (closing rate, divu_adv, opning)
176            !-----------------------------------------------------------------------------!
177            !
178            ! 2.1 closing_net
179            !-----------------
180            ! Compute the net rate of closing due to convergence
181            ! and shear, based on Flato and Hibler (1995).
182            !
183            ! The energy dissipation rate is equal to the net closing rate
184            ! times the ice strength.
185            !
186            ! NOTE: The NET closing rate is equal to the rate that open water
187            !  area is removed, plus the rate at which ice area is removed by
188            !  ridging, minus the rate at which area is added in new ridges.
189            !  The GROSS closing rate is equal to the first two terms (open
190            !  water closing and thin ice ridging) without the third term
191            !  (thick, newly ridged ice).
192
193            closing_net(ji,jj) = rn_cs * 0.5 * ( delta_i(ji,jj) - ABS( divu_i(ji,jj) ) ) - MIN( divu_i(ji,jj), 0._wp )
194
195            ! 2.2 divu_adv
196            !--------------
197            ! Compute divu_adv, the divergence rate given by the transport/
198            ! advection scheme, which may not be equal to divu as computed
199            ! from the velocity field.
200            !
201            ! If divu_adv < 0, make sure the closing rate is large enough
202            ! to give asum = 1.0 after ridging.
203
204            divu_adv(ji,jj) = ( kamax - asum(ji,jj) ) * r1_rdtice  ! asum found in ridgeprep
205
206            IF( divu_adv(ji,jj) < 0._wp )   closing_net(ji,jj) = MAX( closing_net(ji,jj), -divu_adv(ji,jj) )
207
208            ! 2.3 opning
209            !------------
210            ! Compute the (non-negative) opening rate that will give asum = 1.0 after ridging.
211            opning(ji,jj) = closing_net(ji,jj) + divu_adv(ji,jj)
212         END DO
213      END DO
214
215      !-----------------------------------------------------------------------------!
216      ! 3) Ridging iteration
217      !-----------------------------------------------------------------------------!
218      niter           = 1                 ! iteration counter
219      iterate_ridging = 1
220
221      DO WHILE ( iterate_ridging > 0 .AND. niter < nitermax )
222
223         DO jj = 1, jpj
224            DO ji = 1, jpi
225
226               ! 3.2 closing_gross
227               !-----------------------------------------------------------------------------!
228               ! Based on the ITD of ridging and ridged ice, convert the net
229               !  closing rate to a gross closing rate. 
230               ! NOTE: 0 < aksum <= 1
231               closing_gross(ji,jj) = closing_net(ji,jj) / aksum(ji,jj)
232
233               ! correction to closing rate and opening if closing rate is excessive
234               !---------------------------------------------------------------------
235               ! Reduce the closing rate if more than 100% of the open water
236               ! would be removed.  Reduce the opening rate proportionately.
237               IF ( ato_i(ji,jj) > epsi10 .AND. athorn(ji,jj,0) > 0.0 ) THEN
238                  za = athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice
239                  IF ( za > ato_i(ji,jj)) THEN
240                     zfac = ato_i(ji,jj) / za
241                     closing_gross(ji,jj) = closing_gross(ji,jj) * zfac
242                     opning(ji,jj) = opning(ji,jj) * zfac
243                  ENDIF
244               ENDIF
245
246            END DO
247         END DO
248
249         ! correction to closing rate / opening if excessive ice removal
250         !---------------------------------------------------------------
251         ! Reduce the closing rate if more than 100% of any ice category
252         ! would be removed.  Reduce the opening rate proportionately.
253
254         DO jl = 1, jpl
255            DO jj = 1, jpj
256               DO ji = 1, jpi
257                  IF ( a_i(ji,jj,jl) > epsi10 .AND. athorn(ji,jj,jl) > 0._wp )THEN
258                     za = athorn(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice
259                     IF ( za  >  a_i(ji,jj,jl) ) THEN
260                        zfac = a_i(ji,jj,jl) / za
261                        closing_gross(ji,jj) = closing_gross(ji,jj) * zfac
262                        opning       (ji,jj) = opning       (ji,jj) * zfac
263                     ENDIF
264                  ENDIF
265               END DO
266            END DO
267         END DO
268
269         ! 3.3 Redistribute area, volume, and energy.
270         !-----------------------------------------------------------------------------!
271
272         CALL lim_itd_me_ridgeshift( opning, closing_gross, msnow_mlt, esnow_mlt )
273
274         ! 3.4 Compute total area of ice plus open water after ridging.
275         !-----------------------------------------------------------------------------!
276         ! This is in general not equal to one because of divergence during transport
277         asum(:,:) = ato_i(:,:)
278         DO jl = 1, jpl
279            asum(:,:) = asum(:,:) + a_i(:,:,jl)
280         END DO
281
282         ! 3.5 Do we keep on iterating ???
283         !-----------------------------------------------------------------------------!
284         ! Check whether asum = 1.  If not (because the closing and opening
285         ! rates were reduced above), ridge again with new rates.
286
287         iterate_ridging = 0
288
289         DO jj = 1, jpj
290            DO ji = 1, jpi
291               IF (ABS(asum(ji,jj) - kamax ) < epsi10) THEN
292                  closing_net(ji,jj) = 0._wp
293                  opning     (ji,jj) = 0._wp
294               ELSE
295                  iterate_ridging    = 1
296                  divu_adv   (ji,jj) = ( kamax - asum(ji,jj) ) * r1_rdtice
297                  closing_net(ji,jj) = MAX( 0._wp, -divu_adv(ji,jj) )
298                  opning     (ji,jj) = MAX( 0._wp,  divu_adv(ji,jj) )
299               ENDIF
300            END DO
301         END DO
302
303         IF( lk_mpp )   CALL mpp_max( iterate_ridging )
304
305         ! Repeat if necessary.
306         ! NOTE: If strength smoothing is turned on, the ridging must be
307         ! iterated globally because of the boundary update in the
308         ! smoothing.
309
310         niter = niter + 1
311
312         IF( iterate_ridging == 1 ) THEN
313            IF( niter  >  nitermax ) THEN
314               WRITE(numout,*) ' ALERTE : non-converging ridging scheme '
315               WRITE(numout,*) ' niter, iterate_ridging ', niter, iterate_ridging
316            ENDIF
317            CALL lim_itd_me_ridgeprep
318         ENDIF
319
320      END DO !! on the do while over iter
321
322      !-----------------------------------------------------------------------------!
323      ! 4) Ridging diagnostics
324      !-----------------------------------------------------------------------------!
325      ! Convert ridging rate diagnostics to correct units.
326      ! Update fresh water and heat fluxes due to snow melt.
327      DO jj = 1, jpj
328         DO ji = 1, jpi
329
330            dardg1dt(ji,jj) = dardg1dt(ji,jj) * r1_rdtice
331            dardg2dt(ji,jj) = dardg2dt(ji,jj) * r1_rdtice
332            dvirdgdt(ji,jj) = dvirdgdt(ji,jj) * r1_rdtice
333            opening (ji,jj) = opening (ji,jj) * r1_rdtice
334
335            !-----------------------------------------------------------------------------!
336            ! 5) Heat, salt and freshwater fluxes
337            !-----------------------------------------------------------------------------!
338            wfx_snw(ji,jj) = wfx_snw(ji,jj) + msnow_mlt(ji,jj) * r1_rdtice     ! fresh water source for ocean
339            hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + esnow_mlt(ji,jj) * r1_rdtice     ! heat sink for ocean (<0, W.m-2)
340
341         END DO
342      END DO
343
344      ! Check if there is a ridging error
345      IF( lwp ) THEN
346         DO jj = 1, jpj
347            DO ji = 1, jpi
348               IF( ABS( asum(ji,jj) - kamax)  >  epsi10 ) THEN   ! there is a bug
349                  WRITE(numout,*) ' '
350                  WRITE(numout,*) ' ALERTE : Ridging error: total area = ', asum(ji,jj)
351                  WRITE(numout,*) ' limitd_me '
352                  WRITE(numout,*) ' POINT : ', ji, jj
353                  WRITE(numout,*) ' jpl, a_i, athorn '
354                  WRITE(numout,*) 0, ato_i(ji,jj), athorn(ji,jj,0)
355                  DO jl = 1, jpl
356                     WRITE(numout,*) jl, a_i(ji,jj,jl), athorn(ji,jj,jl)
357                  END DO
358               ENDIF
359            END DO
360         END DO
361      END IF
362
363      ! Conservation check
364      IF ( con_i ) THEN
365         CALL lim_column_sum (jpl,   v_i, vt_i_final)
366         fieldid = ' v_i : limitd_me '
367         CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid) 
368      ENDIF
369
370      ! updates
371      CALL lim_var_glo2eqv
372      CALL lim_var_zapsmall
373      CALL lim_var_agg( 1 ) 
374
375      !-----------------------------------------------------------------------------!
376      ! control prints
377      !-----------------------------------------------------------------------------!
378      IF(ln_ctl) THEN
379         CALL prt_ctl_info(' ')
380         CALL prt_ctl_info(' - Cell values : ')
381         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
382         CALL prt_ctl(tab2d_1=e12t , clinfo1=' lim_itd_me  : cell area :')
383         CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_itd_me  : at_i      :')
384         CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_itd_me  : vt_i      :')
385         CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_itd_me  : vt_s      :')
386         DO jl = 1, jpl
387            CALL prt_ctl_info(' ')
388            CALL prt_ctl_info(' - Category : ', ivar1=jl)
389            CALL prt_ctl_info('   ~~~~~~~~~~')
390            CALL prt_ctl(tab2d_1=a_i   (:,:,jl)   , clinfo1= ' lim_itd_me  : a_i      : ')
391            CALL prt_ctl(tab2d_1=ht_i  (:,:,jl)   , clinfo1= ' lim_itd_me  : ht_i     : ')
392            CALL prt_ctl(tab2d_1=ht_s  (:,:,jl)   , clinfo1= ' lim_itd_me  : ht_s     : ')
393            CALL prt_ctl(tab2d_1=v_i   (:,:,jl)   , clinfo1= ' lim_itd_me  : v_i      : ')
394            CALL prt_ctl(tab2d_1=v_s   (:,:,jl)   , clinfo1= ' lim_itd_me  : v_s      : ')
395            CALL prt_ctl(tab2d_1=e_s   (:,:,1,jl) , clinfo1= ' lim_itd_me  : e_s      : ')
396            CALL prt_ctl(tab2d_1=t_su  (:,:,jl)   , clinfo1= ' lim_itd_me  : t_su     : ')
397            CALL prt_ctl(tab2d_1=t_s   (:,:,1,jl) , clinfo1= ' lim_itd_me  : t_snow   : ')
398            CALL prt_ctl(tab2d_1=sm_i  (:,:,jl)   , clinfo1= ' lim_itd_me  : sm_i     : ')
399            CALL prt_ctl(tab2d_1=smv_i (:,:,jl)   , clinfo1= ' lim_itd_me  : smv_i    : ')
400            DO jk = 1, nlay_i
401               CALL prt_ctl_info(' ')
402               CALL prt_ctl_info(' - Layer : ', ivar1=jk)
403               CALL prt_ctl_info('   ~~~~~~~')
404               CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_itd_me  : t_i      : ')
405               CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_itd_me  : e_i      : ')
406            END DO
407         END DO
408      ENDIF
409
410      ! conservation test
411      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limitd_me', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
412
413      ENDIF  ! ln_limdyn=.true.
414      !
415      CALL wrk_dealloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final )
416      !
417      IF( nn_timing == 1 )  CALL timing_stop('limitd_me')
418   END SUBROUTINE lim_itd_me
419
420
421   SUBROUTINE lim_itd_me_icestrength( kstrngth )
422      !!----------------------------------------------------------------------
423      !!                ***  ROUTINE lim_itd_me_icestrength ***
424      !!
425      !! ** Purpose :   computes ice strength used in dynamics routines of ice thickness
426      !!
427      !! ** Method  :   Compute the strength of the ice pack, defined as the energy (J m-2)
428      !!              dissipated per unit area removed from the ice pack under compression,
429      !!              and assumed proportional to the change in potential energy caused
430      !!              by ridging. Note that only Hibler's formulation is stable and that
431      !!              ice strength has to be smoothed
432      !!
433      !! ** Inputs / Ouputs : kstrngth (what kind of ice strength we are using)
434      !!----------------------------------------------------------------------
435      INTEGER, INTENT(in) ::   kstrngth    ! = 1 for Rothrock formulation, 0 for Hibler (1979)
436      INTEGER             ::   ji,jj, jl   ! dummy loop indices
437      INTEGER             ::   ksmooth     ! smoothing the resistance to deformation
438      INTEGER             ::   numts_rm    ! number of time steps for the P smoothing
439      REAL(wp)            ::   zhi, zp, z1_3  ! local scalars
440      REAL(wp), POINTER, DIMENSION(:,:) ::   zworka   ! temporary array used here
441      !!----------------------------------------------------------------------
442
443      CALL wrk_alloc( jpi, jpj, zworka )
444
445      !------------------------------------------------------------------------------!
446      ! 1) Initialize
447      !------------------------------------------------------------------------------!
448      strength(:,:) = 0._wp
449
450      !------------------------------------------------------------------------------!
451      ! 2) Compute thickness distribution of ridging and ridged ice
452      !------------------------------------------------------------------------------!
453      CALL lim_itd_me_ridgeprep
454
455      !------------------------------------------------------------------------------!
456      ! 3) Rothrock(1975)'s method
457      !------------------------------------------------------------------------------!
458      IF( kstrngth == 1 ) THEN
459         z1_3 = 1._wp / 3._wp
460         DO jl = 1, jpl
461            DO jj= 1, jpj
462               DO ji = 1, jpi
463                  !
464                  IF( a_i(ji,jj,jl) > epsi10 .AND. athorn(ji,jj,jl) > 0._wp ) THEN
465                     zhi = v_i(ji,jj,jl) / a_i(ji,jj,jl)
466                     !----------------------------
467                     ! PE loss from deforming ice
468                     !----------------------------
469                     strength(ji,jj) = strength(ji,jj) - athorn(ji,jj,jl) * zhi * zhi
470
471                     !--------------------------
472                     ! PE gain from rafting ice
473                     !--------------------------
474                     strength(ji,jj) = strength(ji,jj) + 2._wp * araft(ji,jj,jl) * zhi * zhi
475
476                     !----------------------------
477                     ! PE gain from ridging ice
478                     !----------------------------
479                     strength(ji,jj) = strength(ji,jj) + aridge(ji,jj,jl) / krdg(ji,jj,jl)     &
480                        * z1_3 * ( hrmax(ji,jj,jl)**2 + hrmin(ji,jj,jl)**2 +  hrmax(ji,jj,jl) * hrmin(ji,jj,jl) ) 
481                        !!(a**3-b**3)/(a-b) = a*a+ab+b*b                     
482                  ENDIF
483                  !
484               END DO
485            END DO
486         END DO
487   
488         strength(:,:) = rn_pe_rdg * Cp * strength(:,:) / aksum(:,:)
489                         ! where Cp = (g/2)*(rhow-rhoi)*(rhoi/rhow) and rn_pe_rdg accounts for frictional dissipation
490         ksmooth = 1
491
492         !------------------------------------------------------------------------------!
493         ! 4) Hibler (1979)' method
494         !------------------------------------------------------------------------------!
495      ELSE                      ! kstrngth ne 1:  Hibler (1979) form
496         !
497         strength(:,:) = rn_pstar * vt_i(:,:) * EXP( - rn_crhg * ( 1._wp - at_i(:,:) )  )
498         !
499         ksmooth = 1
500         !
501      ENDIF                     ! kstrngth
502
503      !
504      !------------------------------------------------------------------------------!
505      ! 5) Impact of brine volume
506      !------------------------------------------------------------------------------!
507      ! CAN BE REMOVED
508      !
509      IF( ln_icestr_bvf ) THEN
510
511         DO jj = 1, jpj
512            DO ji = 1, jpi
513               strength(ji,jj) = strength(ji,jj) * exp(-5.88*SQRT(MAX(bv_i(ji,jj),0.0)))
514            END DO
515         END DO
516
517      ENDIF
518
519      !
520      !------------------------------------------------------------------------------!
521      ! 6) Smoothing ice strength
522      !------------------------------------------------------------------------------!
523      !
524      !-------------------
525      ! Spatial smoothing
526      !-------------------
527      IF ( ksmooth == 1 ) THEN
528
529         CALL lbc_lnk( strength, 'T', 1. )
530
531         DO jj = 2, jpjm1
532            DO ji = 2, jpim1
533               IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > epsi10) THEN ! ice is present
534                  zworka(ji,jj) = ( 4.0 * strength(ji,jj)              &
535                     &                  + strength(ji-1,jj) * tmask(ji-1,jj,1) + strength(ji+1,jj) * tmask(ji+1,jj,1) & 
536                     &                  + strength(ji,jj-1) * tmask(ji,jj-1,1) + strength(ji,jj+1) * tmask(ji,jj+1,1) &
537                     &            ) / ( 4.0 + tmask(ji-1,jj,1) + tmask(ji+1,jj,1) + tmask(ji,jj-1,1) + tmask(ji,jj+1,1) )
538               ELSE
539                  zworka(ji,jj) = 0._wp
540               ENDIF
541            END DO
542         END DO
543
544         DO jj = 2, jpjm1
545            DO ji = 2, jpim1
546               strength(ji,jj) = zworka(ji,jj)
547            END DO
548         END DO
549         CALL lbc_lnk( strength, 'T', 1. )
550
551      ENDIF
552
553      !--------------------
554      ! Temporal smoothing
555      !--------------------
556      IF ( numit == nit000 + nn_fsbc - 1 ) THEN
557         strp1(:,:) = 0.0           
558         strp2(:,:) = 0.0           
559      ENDIF
560
561      IF ( ksmooth == 2 ) THEN
562
563
564         CALL lbc_lnk( strength, 'T', 1. )
565
566         DO jj = 1, jpj - 1
567            DO ji = 1, jpi - 1
568               IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > epsi10) THEN       ! ice is present
569                  numts_rm = 1 ! number of time steps for the running mean
570                  IF ( strp1(ji,jj) > 0.0 ) numts_rm = numts_rm + 1
571                  IF ( strp2(ji,jj) > 0.0 ) numts_rm = numts_rm + 1
572                  zp = ( strength(ji,jj) + strp1(ji,jj) + strp2(ji,jj) ) / numts_rm
573                  strp2(ji,jj) = strp1(ji,jj)
574                  strp1(ji,jj) = strength(ji,jj)
575                  strength(ji,jj) = zp
576
577               ENDIF
578            END DO
579         END DO
580
581      ENDIF ! ksmooth
582
583      CALL lbc_lnk( strength, 'T', 1. )      ! Boundary conditions
584
585      CALL wrk_dealloc( jpi, jpj, zworka )
586      !
587   END SUBROUTINE lim_itd_me_icestrength
588
589
590   SUBROUTINE lim_itd_me_ridgeprep
591      !!---------------------------------------------------------------------!
592      !!                ***  ROUTINE lim_itd_me_ridgeprep ***
593      !!
594      !! ** Purpose :   preparation for ridging and strength calculations
595      !!
596      !! ** Method  :   Compute the thickness distribution of the ice and open water
597      !!              participating in ridging and of the resulting ridges.
598      !!---------------------------------------------------------------------!
599      INTEGER ::   ji,jj, jl    ! dummy loop indices
600      REAL(wp) ::   Gstari, astari, zhi, hrmean, zdummy   ! local scalar
601      REAL(wp), POINTER, DIMENSION(:,:)   ::   zworka    ! temporary array used here
602      REAL(wp), POINTER, DIMENSION(:,:,:) ::   Gsum      ! Gsum(n) = sum of areas in categories 0 to n
603      !------------------------------------------------------------------------------!
604
605      CALL wrk_alloc( jpi,jpj, zworka )
606      CALL wrk_alloc( jpi,jpj,jpl+2, Gsum, kkstart = -1 )
607
608      Gstari     = 1.0/rn_gstar   
609      astari     = 1.0/rn_astar   
610      aksum(:,:)    = 0.0
611      athorn(:,:,:) = 0.0
612      aridge(:,:,:) = 0.0
613      araft (:,:,:) = 0.0
614      hrmin(:,:,:)  = 0.0 
615      hrmax(:,:,:)  = 0.0 
616      hraft(:,:,:)  = 0.0 
617      krdg (:,:,:)  = 1.0
618
619      !     ! Zero out categories with very small areas
620      CALL lim_var_zapsmall
621
622      !------------------------------------------------------------------------------!
623      ! 1) Participation function
624      !------------------------------------------------------------------------------!
625
626      ! Compute total area of ice plus open water.
627      ! This is in general not equal to one because of divergence during transport
628      asum(:,:) = ato_i(:,:)
629      DO jl = 1, jpl
630         asum(:,:) = asum(:,:) + a_i(:,:,jl)
631      END DO
632
633      ! Compute cumulative thickness distribution function
634      ! Compute the cumulative thickness distribution function Gsum,
635      ! where Gsum(n) is the fractional area in categories 0 to n.
636      ! initial value (in h = 0) equals open water area
637
638      Gsum(:,:,-1) = 0._wp
639
640      DO jj = 1, jpj
641         DO ji = 1, jpi
642            IF( ato_i(ji,jj) > epsi10 ) THEN   ;   Gsum(ji,jj,0) = ato_i(ji,jj)
643            ELSE                               ;   Gsum(ji,jj,0) = 0._wp
644            ENDIF
645         END DO
646      END DO
647
648      ! for each value of h, you have to add ice concentration then
649      DO jl = 1, jpl
650         DO jj = 1, jpj 
651            DO ji = 1, jpi
652               IF( a_i(ji,jj,jl) > epsi10 ) THEN   ;   Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1) + a_i(ji,jj,jl)
653               ELSE                                ;   Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1)
654               ENDIF
655            END DO
656         END DO
657      END DO
658
659      ! Normalize the cumulative distribution to 1
660      zworka(:,:) = 1._wp / Gsum(:,:,jpl)
661      DO jl = 0, jpl
662         Gsum(:,:,jl) = Gsum(:,:,jl) * zworka(:,:)
663      END DO
664
665      ! 1.3 Compute participation function a(h) = b(h).g(h) (athorn)
666      !--------------------------------------------------------------------------------------------------
667      ! Compute the participation function athorn; this is analogous to
668      ! a(h) = b(h)g(h) as defined in Thorndike et al. (1975).
669      ! area lost from category n due to ridging/closing
670      ! athorn(n)   = total area lost due to ridging/closing
671      ! assume b(h) = (2/Gstar) * (1 - G(h)/Gstar).
672      !
673      ! The expressions for athorn are found by integrating b(h)g(h) between
674      ! the category boundaries.
675      !-----------------------------------------------------------------
676
677      IF( nn_partfun == 0 ) THEN       !--- Linear formulation (Thorndike et al., 1975)
678         DO jl = 0, jpl   
679            DO jj = 1, jpj 
680               DO ji = 1, jpi
681                  IF( Gsum(ji,jj,jl) < rn_gstar) THEN
682                     athorn(ji,jj,jl) = Gstari * ( Gsum(ji,jj,jl) - Gsum(ji,jj,jl-1) ) * &
683                        &                        ( 2.0 - (Gsum(ji,jj,jl-1) + Gsum(ji,jj,jl) ) * Gstari )
684                  ELSEIF (Gsum(ji,jj,jl-1) < rn_gstar) THEN
685                     athorn(ji,jj,jl) = Gstari * ( rn_gstar - Gsum(ji,jj,jl-1) ) *  &
686                        &                        ( 2.0 - ( Gsum(ji,jj,jl-1) + rn_gstar ) * Gstari )
687                  ELSE
688                     athorn(ji,jj,jl) = 0.0
689                  ENDIF
690               END DO
691            END DO
692         END DO
693
694      ELSE                             !--- Exponential, more stable formulation (Lipscomb et al, 2007)
695         !                       
696         zdummy = 1._wp / ( 1._wp - EXP(-astari) )        ! precompute exponential terms using Gsum as a work array
697         DO jl = -1, jpl
698            Gsum(:,:,jl) = EXP( -Gsum(:,:,jl) * astari ) * zdummy
699         END DO
700         DO jl = 0, jpl
701             athorn(:,:,jl) = Gsum(:,:,jl-1) - Gsum(:,:,jl)
702         END DO
703         !
704      ENDIF
705
706      IF( ln_rafting ) THEN      ! Ridging and rafting ice participation functions
707         !
708         DO jl = 1, jpl
709            DO jj = 1, jpj 
710               DO ji = 1, jpi
711                  IF ( athorn(ji,jj,jl) > 0._wp ) THEN
712!!gm  TANH( -X ) = - TANH( X )  so can be computed only 1 time....
713                     aridge(ji,jj,jl) = ( TANH (  rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) ) + 1.0 ) * 0.5 * athorn(ji,jj,jl)
714                     araft (ji,jj,jl) = ( TANH ( -rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) ) + 1.0 ) * 0.5 * athorn(ji,jj,jl)
715                     IF ( araft(ji,jj,jl) < epsi06 )   araft(ji,jj,jl)  = 0._wp
716                     aridge(ji,jj,jl) = MAX( athorn(ji,jj,jl) - araft(ji,jj,jl), 0.0 )
717                  ENDIF
718               END DO
719            END DO
720         END DO
721
722      ELSE
723         !
724         DO jl = 1, jpl
725            aridge(:,:,jl) = athorn(:,:,jl)
726         END DO
727         !
728      ENDIF
729
730      IF( ln_rafting ) THEN
731
732         IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) > epsi10 .AND. lwp ) THEN
733            DO jl = 1, jpl
734               DO jj = 1, jpj
735                  DO ji = 1, jpi
736                     IF ( aridge(ji,jj,jl) + araft(ji,jj,jl) - athorn(ji,jj,jl) > epsi10 ) THEN
737                        WRITE(numout,*) ' ALERTE 96 : wrong participation function ... '
738                        WRITE(numout,*) ' ji, jj, jl : ', ji, jj, jl
739                        WRITE(numout,*) ' lat, lon   : ', gphit(ji,jj), glamt(ji,jj)
740                        WRITE(numout,*) ' aridge     : ', aridge(ji,jj,1:jpl)
741                        WRITE(numout,*) ' araft      : ', araft(ji,jj,1:jpl)
742                        WRITE(numout,*) ' athorn     : ', athorn(ji,jj,1:jpl)
743                     ENDIF
744                  END DO
745               END DO
746            END DO
747         ENDIF
748
749      ENDIF
750
751      !-----------------------------------------------------------------
752      ! 2) Transfer function
753      !-----------------------------------------------------------------
754      ! Compute max and min ridged ice thickness for each ridging category.
755      ! Assume ridged ice is uniformly distributed between hrmin and hrmax.
756      !
757      ! This parameterization is a modified version of Hibler (1980).
758      ! The mean ridging thickness, hrmean, is proportional to hi^(0.5)
759      !  and for very thick ridging ice must be >= krdgmin*hi
760      !
761      ! The minimum ridging thickness, hrmin, is equal to 2*hi
762      !  (i.e., rafting) and for very thick ridging ice is
763      !  constrained by hrmin <= (hrmean + hi)/2.
764      !
765      ! The maximum ridging thickness, hrmax, is determined by
766      !  hrmean and hrmin.
767      !
768      ! These modifications have the effect of reducing the ice strength
769      ! (relative to the Hibler formulation) when very thick ice is
770      ! ridging.
771      !
772      ! aksum = net area removed/ total area removed
773      ! where total area removed = area of ice that ridges
774      !         net area removed = total area removed - area of new ridges
775      !-----------------------------------------------------------------
776
777      ! Transfer function
778      DO jl = 1, jpl !all categories have a specific transfer function
779         DO jj = 1, jpj
780            DO ji = 1, jpi
781
782               IF (a_i(ji,jj,jl) > epsi10 .AND. athorn(ji,jj,jl) > 0.0 ) THEN
783                  zhi = v_i(ji,jj,jl) / a_i(ji,jj,jl)
784                  hrmean          = MAX(SQRT(rn_hstar*zhi), zhi*krdgmin)
785                  hrmin(ji,jj,jl) = MIN(2.0*zhi, 0.5*(hrmean + zhi))
786                  hrmax(ji,jj,jl) = 2.0*hrmean - hrmin(ji,jj,jl)
787                  hraft(ji,jj,jl) = kraft*zhi
788                  krdg(ji,jj,jl)  = hrmean / zhi
789               ELSE
790                  hraft(ji,jj,jl) = 0.0
791                  hrmin(ji,jj,jl) = 0.0 
792                  hrmax(ji,jj,jl) = 0.0 
793                  krdg (ji,jj,jl) = 1.0
794               ENDIF
795
796            END DO
797         END DO
798      END DO
799
800      ! Normalization factor : aksum, ensures mass conservation
801      aksum(:,:) = athorn(:,:,0)
802      DO jl = 1, jpl 
803         aksum(:,:)    = aksum(:,:) + aridge(:,:,jl) * ( 1._wp - 1._wp / krdg(:,:,jl) )    &
804            &                       + araft (:,:,jl) * ( 1._wp - 1._wp / kraft        )
805      END DO
806      !
807      CALL wrk_dealloc( jpi,jpj, zworka )
808      CALL wrk_dealloc( jpi,jpj,jpl+2, Gsum, kkstart = -1 )
809      !
810   END SUBROUTINE lim_itd_me_ridgeprep
811
812
813   SUBROUTINE lim_itd_me_ridgeshift( opning, closing_gross, msnow_mlt, esnow_mlt )
814      !!----------------------------------------------------------------------
815      !!                ***  ROUTINE lim_itd_me_icestrength ***
816      !!
817      !! ** Purpose :   shift ridging ice among thickness categories of ice thickness
818      !!
819      !! ** Method  :   Remove area, volume, and energy from each ridging category
820      !!              and add to thicker ice categories.
821      !!----------------------------------------------------------------------
822      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   opning         ! rate of opening due to divergence/shear
823      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   closing_gross  ! rate at which area removed, excluding area of new ridges
824      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   msnow_mlt      ! mass of snow added to ocean (kg m-2)
825      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   esnow_mlt      ! energy needed to melt snow in ocean (J m-2)
826      !
827      CHARACTER (len=80) ::   fieldid   ! field identifier
828      LOGICAL, PARAMETER ::   l_conservation_check = .true.  ! if true, check conservation (useful for debugging)
829      !
830      LOGICAL ::   neg_ato_i      ! flag for ato_i(i,j) < -puny
831      LOGICAL ::   large_afrac    ! flag for afrac > 1
832      LOGICAL ::   large_afrft    ! flag for afrac > 1
833      INTEGER ::   ji, jj, jl, jl1, jl2, jk   ! dummy loop indices
834      INTEGER ::   ij                ! horizontal index, combines i and j loops
835      INTEGER ::   icells            ! number of cells with aicen > puny
836      REAL(wp) ::   hL, hR, farea, ztmelts    ! left and right limits of integration
837
838      INTEGER , POINTER, DIMENSION(:) ::   indxi, indxj   ! compressed indices
839
840      REAL(wp), POINTER, DIMENSION(:,:) ::   vice_init, vice_final   ! ice volume summed over categories
841      REAL(wp), POINTER, DIMENSION(:,:) ::   eice_init, eice_final   ! ice energy summed over layers
842
843      REAL(wp), POINTER, DIMENSION(:,:,:) ::   aicen_init, vicen_init   ! ice  area    & volume before ridging
844      REAL(wp), POINTER, DIMENSION(:,:,:) ::   vsnwn_init, esnwn_init   ! snow volume  & energy before ridging
845      REAL(wp), POINTER, DIMENSION(:,:,:) ::   smv_i_init, oa_i_init    ! ice salinity & age    before ridging
846
847      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   eicen_init        ! ice energy before ridging
848
849      REAL(wp), POINTER, DIMENSION(:,:) ::   afrac , fvol     ! fraction of category area ridged & new ridge volume going to n2
850      REAL(wp), POINTER, DIMENSION(:,:) ::   ardg1 , ardg2    ! area of ice ridged & new ridges
851      REAL(wp), POINTER, DIMENSION(:,:) ::   vsrdg , esrdg    ! snow volume & energy of ridging ice
852      REAL(wp), POINTER, DIMENSION(:,:) ::   oirdg1, oirdg2   ! areal age content of ridged & rifging ice
853      REAL(wp), POINTER, DIMENSION(:,:) ::   dhr   , dhr2     ! hrmax - hrmin  &  hrmax^2 - hrmin^2
854
855      REAL(wp), POINTER, DIMENSION(:,:) ::   vrdg1   ! volume of ice ridged
856      REAL(wp), POINTER, DIMENSION(:,:) ::   vrdg2   ! volume of new ridges
857      REAL(wp), POINTER, DIMENSION(:,:) ::   vsw     ! volume of seawater trapped into ridges
858      REAL(wp), POINTER, DIMENSION(:,:) ::   srdg1   ! sal*volume of ice ridged
859      REAL(wp), POINTER, DIMENSION(:,:) ::   srdg2   ! sal*volume of new ridges
860      REAL(wp), POINTER, DIMENSION(:,:) ::   smsw    ! sal*volume of water trapped into ridges
861
862      REAL(wp), POINTER, DIMENSION(:,:) ::   afrft            ! fraction of category area rafted
863      REAL(wp), POINTER, DIMENSION(:,:) ::   arft1 , arft2    ! area of ice rafted and new rafted zone
864      REAL(wp), POINTER, DIMENSION(:,:) ::   virft , vsrft    ! ice & snow volume of rafting ice
865      REAL(wp), POINTER, DIMENSION(:,:) ::   esrft , smrft    ! snow energy & salinity of rafting ice
866      REAL(wp), POINTER, DIMENSION(:,:) ::   oirft1, oirft2   ! areal age content of rafted ice & rafting ice
867
868      REAL(wp), POINTER, DIMENSION(:,:,:) ::   eirft      ! ice energy of rafting ice
869      REAL(wp), POINTER, DIMENSION(:,:,:) ::   erdg1      ! enth*volume of ice ridged
870      REAL(wp), POINTER, DIMENSION(:,:,:) ::   erdg2      ! enth*volume of new ridges
871      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ersw       ! enth of water trapped into ridges
872      !!----------------------------------------------------------------------
873
874      CALL wrk_alloc( (jpi+1)*(jpj+1),      indxi, indxj )
875      CALL wrk_alloc( jpi, jpj,             vice_init, vice_final, eice_init, eice_final )
876      CALL wrk_alloc( jpi, jpj,             afrac, fvol , ardg1, ardg2, vsrdg, esrdg, oirdg1, oirdg2, dhr, dhr2 )
877      CALL wrk_alloc( jpi, jpj,             vrdg1, vrdg2, vsw  , srdg1, srdg2, smsw )
878      CALL wrk_alloc( jpi, jpj,             afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 )
879      CALL wrk_alloc( jpi, jpj, jpl,        aicen_init, vicen_init, vsnwn_init, esnwn_init, smv_i_init, oa_i_init )
880      CALL wrk_alloc( jpi, jpj, nlay_i+1,      eirft, erdg1, erdg2, ersw )
881      CALL wrk_alloc( jpi, jpj, nlay_i+1, jpl, eicen_init )
882
883      ! Conservation check
884      eice_init(:,:) = 0._wp
885
886      IF( con_i ) THEN
887         CALL lim_column_sum        (jpl,    v_i,       vice_init )
888         CALL lim_column_sum_energy (jpl, nlay_i,  e_i, eice_init )
889         DO ji = mi0(iiceprt), mi1(iiceprt)
890            DO jj = mj0(jiceprt), mj1(jiceprt)
891               WRITE(numout,*) ' vice_init  : ', vice_init(ji,jj)
892               WRITE(numout,*) ' eice_init  : ', eice_init(ji,jj)
893            END DO
894         END DO
895      ENDIF
896
897      !-------------------------------------------------------------------------------
898      ! 1) Compute change in open water area due to closing and opening.
899      !-------------------------------------------------------------------------------
900
901      neg_ato_i = .false.
902
903      DO jj = 1, jpj
904         DO ji = 1, jpi
905            ato_i(ji,jj) = ato_i(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice        &
906               &                        + opning(ji,jj)                          * rdt_ice
907            IF( ato_i(ji,jj) < -epsi10 ) THEN
908               neg_ato_i = .TRUE.
909            ELSEIF( ato_i(ji,jj) < 0._wp ) THEN    ! roundoff error
910               ato_i(ji,jj) = 0._wp
911            ENDIF
912         END DO
913      END DO
914
915      ! if negative open water area alert it
916      IF( neg_ato_i .AND. lwp ) THEN       ! there is a bug
917         DO jj = 1, jpj 
918            DO ji = 1, jpi
919               IF( ato_i(ji,jj) < -epsi10 ) THEN
920                  WRITE(numout,*) '' 
921                  WRITE(numout,*) 'Ridging error: ato_i < 0'
922                  WRITE(numout,*) 'ato_i : ', ato_i(ji,jj)
923               ENDIF
924            END DO
925         END DO
926      ENDIF
927
928      !-----------------------------------------------------------------
929      ! 2) Save initial state variables
930      !-----------------------------------------------------------------
931
932      DO jl = 1, jpl
933         aicen_init(:,:,jl) = a_i(:,:,jl)
934         vicen_init(:,:,jl) = v_i(:,:,jl)
935         vsnwn_init(:,:,jl) = v_s(:,:,jl)
936         !
937         smv_i_init(:,:,jl) = smv_i(:,:,jl)
938         oa_i_init (:,:,jl) = oa_i (:,:,jl)
939      END DO
940
941      esnwn_init(:,:,:) = e_s(:,:,1,:)
942
943      DO jl = 1, jpl 
944         DO jk = 1, nlay_i
945            eicen_init(:,:,jk,jl) = e_i(:,:,jk,jl)
946         END DO
947      END DO
948
949      !
950      !-----------------------------------------------------------------
951      ! 3) Pump everything from ice which is being ridged / rafted
952      !-----------------------------------------------------------------
953      ! Compute the area, volume, and energy of ice ridging in each
954      ! category, along with the area of the resulting ridge.
955
956      DO jl1 = 1, jpl !jl1 describes the ridging category
957
958         !------------------------------------------------
959         ! 3.1) Identify grid cells with nonzero ridging
960         !------------------------------------------------
961
962         icells = 0
963         DO jj = 1, jpj
964            DO ji = 1, jpi
965               IF( aicen_init(ji,jj,jl1) > epsi10 .AND. athorn(ji,jj,jl1) > 0._wp  &
966                  &   .AND. closing_gross(ji,jj) > 0._wp ) THEN
967                  icells = icells + 1
968                  indxi(icells) = ji
969                  indxj(icells) = jj
970               ENDIF
971            END DO
972         END DO
973
974         large_afrac = .false.
975         large_afrft = .false.
976
977         DO ij = 1, icells
978            ji = indxi(ij)
979            jj = indxj(ij)
980
981            !--------------------------------------------------------------------
982            ! 3.2) Compute area of ridging ice (ardg1) and of new ridge (ardg2)
983            !--------------------------------------------------------------------
984
985            ardg1(ji,jj) = aridge(ji,jj,jl1)*closing_gross(ji,jj)*rdt_ice
986            arft1(ji,jj) = araft (ji,jj,jl1)*closing_gross(ji,jj)*rdt_ice
987            ardg2(ji,jj) = ardg1(ji,jj) / krdg(ji,jj,jl1)
988            arft2(ji,jj) = arft1(ji,jj) / kraft
989
990            oirdg1(ji,jj)= aridge(ji,jj,jl1)*closing_gross(ji,jj)*rdt_ice
991            oirft1(ji,jj)= araft (ji,jj,jl1)*closing_gross(ji,jj)*rdt_ice
992            oirdg2(ji,jj)= oirdg1(ji,jj) / krdg(ji,jj,jl1)
993            oirft2(ji,jj)= oirft1(ji,jj) / kraft
994
995            !---------------------------------------------------------------
996            ! 3.3) Compute ridging /rafting fractions, make sure afrac <=1
997            !---------------------------------------------------------------
998
999            afrac(ji,jj) = ardg1(ji,jj) / aicen_init(ji,jj,jl1) !ridging
1000            afrft(ji,jj) = arft1(ji,jj) / aicen_init(ji,jj,jl1) !rafting
1001
1002            IF (afrac(ji,jj) > kamax + epsi10) THEN  !riging
1003               large_afrac = .true.
1004            ELSEIF (afrac(ji,jj) > kamax) THEN  ! roundoff error
1005               afrac(ji,jj) = kamax
1006            ENDIF
1007            IF (afrft(ji,jj) > kamax + epsi10) THEN !rafting
1008               large_afrft = .true.
1009            ELSEIF (afrft(ji,jj) > kamax) THEN  ! roundoff error
1010               afrft(ji,jj) = kamax
1011            ENDIF
1012
1013            !--------------------------------------------------------------------------
1014            ! 3.4) Subtract area, volume, and energy from ridging
1015            !     / rafting category n1.
1016            !--------------------------------------------------------------------------
1017            vrdg1(ji,jj) = vicen_init(ji,jj,jl1) * afrac(ji,jj)
1018            vrdg2(ji,jj) = vrdg1(ji,jj) * ( 1. + rn_por_rdg )
1019            vsw  (ji,jj) = vrdg1(ji,jj) * rn_por_rdg
1020
1021            vsrdg(ji,jj) = vsnwn_init(ji,jj,jl1) * afrac(ji,jj)
1022            esrdg(ji,jj) = esnwn_init(ji,jj,jl1) * afrac(ji,jj)
1023            srdg1(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj)
1024            srdg2(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) !! MV HC 2014 this line seems useless
1025
1026            ! rafting volumes, heat contents ...
1027            virft(ji,jj) = vicen_init(ji,jj,jl1) * afrft(ji,jj)
1028            vsrft(ji,jj) = vsnwn_init(ji,jj,jl1) * afrft(ji,jj)
1029            esrft(ji,jj) = esnwn_init(ji,jj,jl1) * afrft(ji,jj)
1030            smrft(ji,jj) = smv_i_init(ji,jj,jl1) * afrft(ji,jj) 
1031
1032            ! substract everything
1033            a_i(ji,jj,jl1)   = a_i(ji,jj,jl1)   - ardg1(ji,jj)  - arft1(ji,jj)
1034            v_i(ji,jj,jl1)   = v_i(ji,jj,jl1)   - vrdg1(ji,jj)  - virft(ji,jj)
1035            v_s(ji,jj,jl1)   = v_s(ji,jj,jl1)   - vsrdg(ji,jj)  - vsrft(ji,jj)
1036            e_s(ji,jj,1,jl1) = e_s(ji,jj,1,jl1) - esrdg(ji,jj)  - esrft(ji,jj)
1037            oa_i(ji,jj,jl1)  = oa_i(ji,jj,jl1)  - oirdg1(ji,jj) - oirft1(ji,jj)
1038            smv_i(ji,jj,jl1) = smv_i(ji,jj,jl1) - srdg1(ji,jj)  - smrft(ji,jj)
1039
1040            !-----------------------------------------------------------------
1041            ! 3.5) Compute properties of new ridges
1042            !-----------------------------------------------------------------
1043            !-------------
1044            ! Salinity
1045            !-------------
1046            smsw(ji,jj)  = vsw(ji,jj) * sss_m(ji,jj)                      ! salt content of seawater frozen in voids !! MV HC2014
1047            srdg2(ji,jj) = srdg1(ji,jj) + smsw(ji,jj)                     ! salt content of new ridge
1048
1049            !srdg2(ji,jj) = MIN( rn_simax * vrdg2(ji,jj) , zsrdg2 )         ! impose a maximum salinity
1050           
1051            sfx_dyn(ji,jj) = sfx_dyn(ji,jj) - smsw(ji,jj) * rhoic * r1_rdtice
1052            wfx_dyn(ji,jj) = wfx_dyn(ji,jj) - vsw (ji,jj) * rhoic * r1_rdtice   ! gurvan: increase in ice volume du to seawater frozen in voids             
1053
1054            !------------------------------------           
1055            ! 3.6 Increment ridging diagnostics
1056            !------------------------------------           
1057
1058            !        jl1 looping 1-jpl
1059            !           ij looping 1-icells
1060
1061            dardg1dt   (ji,jj) = dardg1dt(ji,jj) + ardg1(ji,jj) + arft1(ji,jj)
1062            dardg2dt   (ji,jj) = dardg2dt(ji,jj) + ardg2(ji,jj) + arft2(ji,jj)
1063            opening    (ji,jj) = opening (ji,jj) + opning(ji,jj) * rdt_ice
1064
1065            IF( con_i )   vice_init(ji,jj) = vice_init(ji,jj) + vrdg2(ji,jj) - vrdg1(ji,jj)
1066
1067            !------------------------------------------           
1068            ! 3.7 Put the snow somewhere in the ocean
1069            !------------------------------------------           
1070            !  Place part of the snow lost by ridging into the ocean.
1071            !  Note that esnow_mlt < 0; the ocean must cool to melt snow.
1072            !  If the ocean temp = Tf already, new ice must grow.
1073            !  During the next time step, thermo_rates will determine whether
1074            !  the ocean cools or new ice grows.
1075            !        jl1 looping 1-jpl
1076            !           ij looping 1-icells
1077
1078            msnow_mlt(ji,jj) = msnow_mlt(ji,jj) + rhosn*vsrdg(ji,jj)*(1.0-rn_fsnowrdg)   &   ! rafting included
1079               &                                + rhosn*vsrft(ji,jj)*(1.0-rn_fsnowrft)
1080
1081            ! in J/m2 (same as e_s)
1082            esnow_mlt(ji,jj) = esnow_mlt(ji,jj) - esrdg(ji,jj)*(1.0-rn_fsnowrdg)         &   !rafting included
1083               &                                - esrft(ji,jj)*(1.0-rn_fsnowrft)         
1084
1085            !-----------------------------------------------------------------
1086            ! 3.8 Compute quantities used to apportion ice among categories
1087            ! in the n2 loop below
1088            !-----------------------------------------------------------------
1089
1090            !        jl1 looping 1-jpl
1091            !           ij looping 1-icells
1092
1093            dhr (ji,jj) = hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1)
1094            dhr2(ji,jj) = hrmax(ji,jj,jl1) * hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) * hrmin(ji,jj,jl1)
1095
1096         END DO
1097
1098         !--------------------------------------------------------------------
1099         ! 3.9 Compute ridging ice enthalpy, remove it from ridging ice and
1100         !      compute ridged ice enthalpy
1101         !--------------------------------------------------------------------
1102         DO jk = 1, nlay_i
1103            DO ij = 1, icells
1104               ji = indxi(ij)
1105               jj = indxj(ij)
1106               ! heat content of ridged ice
1107               erdg1(ji,jj,jk)      = eicen_init(ji,jj,jk,jl1) * afrac(ji,jj) 
1108               eirft(ji,jj,jk)      = eicen_init(ji,jj,jk,jl1) * afrft(ji,jj)
1109               e_i  (ji,jj,jk,jl1)  = e_i(ji,jj,jk,jl1) - erdg1(ji,jj,jk) - eirft(ji,jj,jk)
1110               
1111               
1112               ! enthalpy of the trapped seawater (J/m2, >0)
1113               ! clem: if sst>0, then ersw <0 (is that possible?)
1114               ersw(ji,jj,jk)   = - rhoic * vsw(ji,jj) * rcp * sst_m(ji,jj) * r1_nlay_i
1115
1116               ! heat flux to the ocean
1117               hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ersw(ji,jj,jk) * r1_rdtice  ! > 0 [W.m-2] ocean->ice flux
1118
1119               ! it is added to sea ice because the sign convention is the opposite of the sign convention for the ocean
1120               erdg2(ji,jj,jk)  = erdg1(ji,jj,jk) + ersw(ji,jj,jk)
1121
1122            END DO
1123         END DO
1124
1125
1126         IF( con_i ) THEN
1127            DO jk = 1, nlay_i
1128               DO ij = 1, icells
1129                  ji = indxi(ij)
1130                  jj = indxj(ij)
1131                  eice_init(ji,jj) = eice_init(ji,jj) + erdg2(ji,jj,jk) - erdg1(ji,jj,jk)
1132               END DO
1133            END DO
1134         ENDIF
1135
1136         IF( large_afrac .AND. lwp ) THEN   ! there is a bug
1137            DO ij = 1, icells
1138               ji = indxi(ij)
1139               jj = indxj(ij)
1140               IF( afrac(ji,jj) > kamax + epsi10 ) THEN
1141                  WRITE(numout,*) ''
1142                  WRITE(numout,*) ' ardg > a_i'
1143                  WRITE(numout,*) ' ardg, aicen_init : ', ardg1(ji,jj), aicen_init(ji,jj,jl1)
1144               ENDIF
1145            END DO
1146         ENDIF
1147         IF( large_afrft .AND. lwp ) THEN  ! there is a bug
1148            DO ij = 1, icells
1149               ji = indxi(ij)
1150               jj = indxj(ij)
1151               IF( afrft(ji,jj) > kamax + epsi10 ) THEN
1152                  WRITE(numout,*) ''
1153                  WRITE(numout,*) ' arft > a_i'
1154                  WRITE(numout,*) ' arft, aicen_init : ', arft1(ji,jj), aicen_init(ji,jj,jl1)
1155               ENDIF
1156            END DO
1157         ENDIF
1158
1159         !-------------------------------------------------------------------------------
1160         ! 4) Add area, volume, and energy of new ridge to each category jl2
1161         !-------------------------------------------------------------------------------
1162         !        jl1 looping 1-jpl
1163         DO jl2  = 1, jpl 
1164            ! over categories to which ridged ice is transferred
1165            DO ij = 1, icells
1166               ji = indxi(ij)
1167               jj = indxj(ij)
1168
1169               ! Compute the fraction of ridged ice area and volume going to
1170               ! thickness category jl2.
1171               ! Transfer area, volume, and energy accordingly.
1172
1173               IF( hrmin(ji,jj,jl1) >= hi_max(jl2) .OR. hrmax(ji,jj,jl1) <= hi_max(jl2-1) ) THEN
1174                  hL = 0._wp
1175                  hR = 0._wp
1176               ELSE
1177                  hL = MAX( hrmin(ji,jj,jl1), hi_max(jl2-1) )
1178                  hR = MIN( hrmax(ji,jj,jl1), hi_max(jl2)   )
1179               ENDIF
1180
1181               ! fraction of ridged ice area and volume going to n2
1182               farea = ( hR - hL ) / dhr(ji,jj) 
1183               fvol(ji,jj) = ( hR*hR - hL*hL ) / dhr2(ji,jj)
1184
1185               a_i  (ji,jj  ,jl2) = a_i  (ji,jj  ,jl2) + ardg2 (ji,jj) * farea
1186               v_i  (ji,jj  ,jl2) = v_i  (ji,jj  ,jl2) + vrdg2 (ji,jj) * fvol(ji,jj)
1187               v_s  (ji,jj  ,jl2) = v_s  (ji,jj  ,jl2) + vsrdg (ji,jj) * fvol(ji,jj) * rn_fsnowrdg
1188               e_s  (ji,jj,1,jl2) = e_s  (ji,jj,1,jl2) + esrdg (ji,jj) * fvol(ji,jj) * rn_fsnowrdg
1189               smv_i(ji,jj  ,jl2) = smv_i(ji,jj  ,jl2) + srdg2 (ji,jj) * fvol(ji,jj)
1190               oa_i (ji,jj  ,jl2) = oa_i (ji,jj  ,jl2) + oirdg2(ji,jj) * farea
1191
1192            END DO ! ij
1193
1194            ! Transfer ice energy to category jl2 by ridging
1195            DO jk = 1, nlay_i
1196               DO ij = 1, icells
1197                  ji = indxi(ij)
1198                  jj = indxj(ij)
1199                  e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + fvol(ji,jj) * erdg2(ji,jj,jk)
1200               END DO
1201            END DO
1202            !
1203         END DO                 ! jl2 (new ridges)           
1204
1205         DO jl2 = 1, jpl 
1206
1207            DO ij = 1, icells
1208               ji = indxi(ij)
1209               jj = indxj(ij)
1210               ! Compute the fraction of rafted ice area and volume going to
1211               ! thickness category jl2, transfer area, volume, and energy accordingly.
1212               !
1213               IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. hraft(ji,jj,jl1) >  hi_max(jl2-1) ) THEN
1214                  a_i  (ji,jj  ,jl2) = a_i  (ji,jj  ,jl2) + arft2 (ji,jj)
1215                  v_i  (ji,jj  ,jl2) = v_i  (ji,jj  ,jl2) + virft (ji,jj)
1216                  v_s  (ji,jj  ,jl2) = v_s  (ji,jj  ,jl2) + vsrft (ji,jj) * rn_fsnowrft
1217                  e_s  (ji,jj,1,jl2) = e_s  (ji,jj,1,jl2) + esrft (ji,jj) * rn_fsnowrft
1218                  smv_i(ji,jj  ,jl2) = smv_i(ji,jj  ,jl2) + smrft (ji,jj)   
1219                  oa_i (ji,jj  ,jl2) = oa_i (ji,jj  ,jl2) + oirft2(ji,jj)   
1220               ENDIF
1221               !
1222            END DO
1223
1224            ! Transfer rafted ice energy to category jl2
1225            DO jk = 1, nlay_i
1226               DO ij = 1, icells
1227                  ji = indxi(ij)
1228                  jj = indxj(ij)
1229                  IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. hraft(ji,jj,jl1) > hi_max(jl2-1)  ) THEN
1230                     e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + eirft(ji,jj,jk)
1231                  ENDIF
1232               END DO
1233            END DO
1234
1235         END DO
1236
1237      END DO ! jl1 (deforming categories)
1238
1239      ! Conservation check
1240      IF ( con_i ) THEN
1241         CALL lim_column_sum (jpl,   v_i, vice_final)
1242         fieldid = ' v_i : limitd_me '
1243         CALL lim_cons_check (vice_init, vice_final, 1.0e-6, fieldid) 
1244
1245         CALL lim_column_sum_energy (jpl, nlay_i,  e_i, eice_final )
1246         fieldid = ' e_i : limitd_me '
1247         CALL lim_cons_check (eice_init, eice_final, 1.0e-2, fieldid) 
1248
1249         DO ji = mi0(iiceprt), mi1(iiceprt)
1250            DO jj = mj0(jiceprt), mj1(jiceprt)
1251               WRITE(numout,*) ' vice_init  : ', vice_init (ji,jj)
1252               WRITE(numout,*) ' vice_final : ', vice_final(ji,jj)
1253               WRITE(numout,*) ' eice_init  : ', eice_init (ji,jj)
1254               WRITE(numout,*) ' eice_final : ', eice_final(ji,jj)
1255            END DO
1256         END DO
1257      ENDIF
1258      !
1259      CALL wrk_dealloc( (jpi+1)*(jpj+1),      indxi, indxj )
1260      CALL wrk_dealloc( jpi, jpj,             vice_init, vice_final, eice_init, eice_final )
1261      CALL wrk_dealloc( jpi, jpj,             afrac, fvol , ardg1, ardg2, vsrdg, esrdg, oirdg1, oirdg2, dhr, dhr2 )
1262      CALL wrk_dealloc( jpi, jpj,             vrdg1, vrdg2, vsw  , srdg1, srdg2, smsw )
1263      CALL wrk_dealloc( jpi, jpj,             afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 )
1264      CALL wrk_dealloc( jpi, jpj, jpl,        aicen_init, vicen_init, vsnwn_init, esnwn_init, smv_i_init, oa_i_init )
1265      CALL wrk_dealloc( jpi, jpj, nlay_i+1,      eirft, erdg1, erdg2, ersw )
1266      CALL wrk_dealloc( jpi, jpj, nlay_i+1, jpl, eicen_init )
1267      !
1268   END SUBROUTINE lim_itd_me_ridgeshift
1269
1270   SUBROUTINE lim_itd_me_init
1271      !!-------------------------------------------------------------------
1272      !!                   ***  ROUTINE lim_itd_me_init ***
1273      !!
1274      !! ** Purpose :   Physical constants and parameters linked
1275      !!                to the mechanical ice redistribution
1276      !!
1277      !! ** Method  :   Read the namiceitdme namelist
1278      !!                and check the parameters values
1279      !!                called at the first timestep (nit000)
1280      !!
1281      !! ** input   :   Namelist namiceitdme
1282      !!-------------------------------------------------------------------
1283      INTEGER :: ios                 ! Local integer output status for namelist read
1284      NAMELIST/namiceitdme/ rn_cs, rn_fsnowrdg, rn_fsnowrft,              & 
1285        &                   rn_gstar, rn_astar, rn_hstar, ln_rafting, rn_hraft, rn_craft, rn_por_rdg, &
1286        &                   nn_partfun
1287      !!-------------------------------------------------------------------
1288      !
1289      REWIND( numnam_ice_ref )              ! Namelist namicetdme in reference namelist : Ice mechanical ice redistribution
1290      READ  ( numnam_ice_ref, namiceitdme, IOSTAT = ios, ERR = 901)
1291901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitdme in reference namelist', lwp )
1292
1293      REWIND( numnam_ice_cfg )              ! Namelist namiceitdme in configuration namelist : Ice mechanical ice redistribution
1294      READ  ( numnam_ice_cfg, namiceitdme, IOSTAT = ios, ERR = 902 )
1295902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitdme in configuration namelist', lwp )
1296      IF(lwm) WRITE ( numoni, namiceitdme )
1297      !
1298      IF (lwp) THEN                          ! control print
1299         WRITE(numout,*)
1300         WRITE(numout,*)' lim_itd_me_init : ice parameters for mechanical ice redistribution '
1301         WRITE(numout,*)' ~~~~~~~~~~~~~~~'
1302         WRITE(numout,*)'   Fraction of shear energy contributing to ridging        rn_cs       = ', rn_cs 
1303         WRITE(numout,*)'   Fraction of snow volume conserved during ridging        rn_fsnowrdg = ', rn_fsnowrdg 
1304         WRITE(numout,*)'   Fraction of snow volume conserved during ridging        rn_fsnowrft = ', rn_fsnowrft 
1305         WRITE(numout,*)'   Fraction of total ice coverage contributing to ridging  rn_gstar    = ', rn_gstar
1306         WRITE(numout,*)'   Equivalent to G* for an exponential part function       rn_astar    = ', rn_astar
1307         WRITE(numout,*)'   Quantity playing a role in max ridged ice thickness     rn_hstar    = ', rn_hstar
1308         WRITE(numout,*)'   Rafting of ice sheets or not                            ln_rafting  = ', ln_rafting
1309         WRITE(numout,*)'   Parmeter thickness (threshold between ridge-raft)       rn_hraft    = ', rn_hraft
1310         WRITE(numout,*)'   Rafting hyperbolic tangent coefficient                  rn_craft    = ', rn_craft 
1311         WRITE(numout,*)'   Initial porosity of ridges                              rn_por_rdg  = ', rn_por_rdg
1312         WRITE(numout,*)'   Switch for part. function (0) linear (1) exponential    nn_partfun  = ', nn_partfun
1313      ENDIF
1314      !
1315   END SUBROUTINE lim_itd_me_init
1316
1317#else
1318   !!----------------------------------------------------------------------
1319   !!   Default option         Empty module          NO LIM-3 sea-ice model
1320   !!----------------------------------------------------------------------
1321CONTAINS
1322   SUBROUTINE lim_itd_me           ! Empty routines
1323   END SUBROUTINE lim_itd_me
1324   SUBROUTINE lim_itd_me_icestrength
1325   END SUBROUTINE lim_itd_me_icestrength
1326   SUBROUTINE lim_itd_me_init
1327   END SUBROUTINE lim_itd_me_init
1328#endif
1329   !!======================================================================
1330END MODULE limitd_me
Note: See TracBrowser for help on using the repository browser.