source: branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 @ 5086

Last change on this file since 5086 was 5086, checked in by timgraham, 7 years ago

Merged head of trunk into branch in preparation for putting code back onto the trunk
In working copy ran the command:
svn merge svn+sshtimgraham@…/ipsl/forge/projets/nemo/svn/trunk

Also recompiled NEMO_book.pdf with merged input files

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