New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
limitd_me.F90 in trunk/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

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

Last change on this file since 3558 was 3558, checked in by rblod, 11 years ago

Fix issues when using key_nosignedzeo, see ticket #996

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