source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 @ 4298

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